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Summary. The linearized potential equation for unsteady motion in frictionless, 
supersonic flow is transformed from the classical wave equation to the canonical form 
bez — byy — Oss = %-, With the aid of a modified Lorentz transformation. Possible in- 
variant transformations of the latter, including the classical Lorentz transformation, 
are discussed. Eleven coordinate systems (each of which has its counterpart in the classical 
theory of the wave equation) permitting separation of variables are set forth, their 
derivation being based on the analogy between the hyperbolic metric defined by 
(ds)” = (dx)* — (dy)’ — (dz)? and the Euclidean (Cartesian) metric. A few practical 
applications are indicated. 

1. Introduction. We consider here the linearized equation for the velocity potential 
in unsteady, supersonic flow and various coordinate transformations to which it may be 
usefully subjected in order to effect separation of variables. 

In connection with the linearization of the original equations we remark that the 
assumptions implicit therein impose much stronger restrictions than in such classical 
fields as electricity and magnetism. Let 6 be the fineness ratio, defined as the larger of 
the maximum thickness of a wing, or the amplitude of transverse motion, divided by 
the maximum wing chord I, the latter serving as the characteristic length throughout 
the following analysis. Further let » be a dimensionless measure of time rate of change, 
where the characteristic time is lc~*, c being the sonic velocity in the undisturbed medium, 
and let sl be the average wing span. Then an extension of the two dimensional analysis 
of Lin, Reissner and Tsien’ shows that sufficient conditions for linearization (M > 1) are 


Ms <1, vMi<K1 (1.1) 


and any one or more of 
| M—1|,», oc € ef”, (1.2) 


where is the free stream Mach number. 
In the case of slender bodies slightly different restrictions obtain.” 


*Received Apri] 1, 1953. The material in this paper is drawn from lectures given while the author 
was on sabbatical leave in London, particularly at the Imperial College of Science and Technology 
(February and March, 1952). 

IC, C. Lin, E. Reissner and H. 8. Tsien, J. Math. & Ph. 27, 220 (1948). 

2J. W. Miles, J. Aero. Sci. 19, 380 (1952). 
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Because of these restrictions the results obtained herein probably have but limited 
practical application. Nevertheless, we feel that they are of interest per se and, in addition, 
may lead to useful results in the hands of other investigators. Moreover, they may be 
of some value in attacking the non-linear equations. 

2. The potential equation. The linearized potential equation governing small dis- 
turbances with respect to a fixed coordinate system in a perfect fluid is (the wave 
equation) 

oxx + dyy + Ozz = Orr, (2.1) 
where X, Y, Z are dimensionless Cartesian coordinates in a fixed reference frame, T 
is a dimensionless time obtained by multiplying the true time (¢) by the sonic velocity 
(c) and dividing by /, and ¢ is the dimensionless velocity potential (reference quantity: 
Ul). In the case of a supersonic flight at Mach number M along the negative X axis, 
the body in question may be brought to rest and (1) reduced to (what may be regarded 
as) a canonical form by introducing the modified Lorentz transformation 


(2) aaa ( wy) a. ae (2.2) 


in which case we have 
Old = des — Dyy — Des = Der (2.3) 


where 1) may be designated as the “hyperbolic Laplacian” operator.* 

Equation (3) may be identified as the wave equation in the coordinates (x, iy, 7z, 7), 
just as Bateman® has identified the classical wave equation as Laplace’s equation, 
albeit four dimensional, in (X, Y, Z, 77). Such an identification suggests, by analogy, 
numerous solutions to (3), which, to be sure, must be appropriately restricted if they 
are to correspond to physical reality. 

If in (3) we pose the harmonic dependencef exp (—cxr), with the generality implied 


by Fourier’s theorem, we obtain 


D¢+ «¢=0 (2.4) 


which may be appropriately designated as the “hyperbolic Helmholtz equation”’. 

The origin of the moving (xz, y, z) coordinates is conveniently chosen at the most 
upstream point of the body creating the disturbance. Then, in consequence of the 
supersonic flight velocity, we may assert 


@ = 0, we @. (2.5) 
Accordingly, it is often expedient to introduce the Laplace transformation 
@= Lid} = | e “ddx (2.6) 


*This notation has been used by P. A. Lagerstrom in unpublished lectures at the California Institute 
of Technology. 

3H. Bateman, Partial differential equations, Dover Publ., New York, 1944, p. 384. 

+The choice x = (M? — 1)-1/2(@/c) yields the dependence exp (iwt) in the original time domain; 


cf. (2.9) infra. 
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which, applied to (4), yields 
6,, + 6, — 6 = 0, (2.7) 
» = 3° a Fo (2.8) 


In the application of the end results, it is most convenient to deal with a coordinate 
z* measured from the foremost point on the airfoil and the true time ¢, the corresponding 
potential being given by 


o*(z*, ¥, 2, EF) = AZ, ¥, 2, 7) 
= ¢|(M* — 1)7'2x*, y, 2, M(M? — 1)°'72* — (M’ — 1)'°Ct/D). ~~ (2.9) 


However, in all of the subsequent discussion we shall deal implicitly with x and r. 

3. Invariant transformations. The general question of transformations under which 
the classical wave equation remains invariant has been discussed in a series of papers 
by Bateman.* The corresponding transformations of (2.3) and (2.4) follow via the 
analogy suggested above. 

It is immediately evident that (2.3) is invariant under independent translations of 
(x, y, 2, 7), a spherical rotation of (y, z, 7) with x fixed in direction (N.B.: (2.3) is not 
invariant under a rotation involving 2, as, e.g., in the case of a transformation to Mach 
coordinates.), a (simultaneous) scale transformation of (x, y, 2, 7), and a scale trans- 
formation of ¢. 

Rather less obvious are the inversions studied by Bateman. The simplest of these 
imply that if ¢(z, y, z, 7) is a solution to (2.3) so also are 


V(x, y, 2, 7) = ud(ux, wy, wz, wT), (3.1) 
p=(r —-y-—-2- 7)", (3.2) 
x(t, y, 2, 7) = lly +24 7°42), do +24 7 — 2D), v2, v7], (3.3) 
y=(x-—y)'. (3.4) 


As an example of a more general result, a homogeneous solution of degree —1 is given by 


d(x, y,z, 7) = (x ty) (a — y)* (2 + ir) *e - ir)? 


a be 
r—y—2—ry(-l "Pia BY 2 (3.5) 
La’ B’ 7’ 


from which additional transformations may be obtained via the many transformations 
of the generalized hypergeometric function P (in the Riemann-Papperitz notation). 
Perhaps the most interesting (but not necessarily the most important, since many 
valuable inferences are afforded by the various scale transformations) transformation 
under which (2.3) remains invariant is that of Lorentz, which is most conveniently 
written in the normalized form (so that all transformations obtained by assigning 


‘H. Bateman, Proc. London Math. Soc. (2) 8, 223 (1909); ibid 7, 70 (1909); 8, 469 (1910); 10, 
7 (1911). 
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. e . 6 
different values to m are members of a group; cf. the work of Lagerstrom® and Hayes 


in steady flow). 


(7 ] = (1 — m’)" sh m\(2), m|< i. 3.6) 
y m | y 


This result has its most direct application in extending a solution for a wing leaving the 


wing tip y = 0 to one making an angle tan * m with the direction of flight. (But if m 
is negative the original tip becomes a trailing edge, and the Kutta condition then must 
be introduced.) This result has been applied to the problem of a rated wing tip in unsteady 
flow, starting from the known solution for a rectangular wing and applying the Kutta 
condition where the edge is trailing (see J. W. Miles, Q. Appl. Math., 11, 363; 1953). 
Additional transformations under which (2.3) remains invariant may be obtained 
from symmetry considerations among (y, z, 7). Finally, (2.4) is invariant under all of 
the foregoing transformations not involving 7 and to a simultaneous scale transformation 


of (x, y, z, K~'). 


4. Coordinate transformations. We shall consider only those coordinate transforma- 
tions 


L,Y, 2 = L,Y, AM » V2 » Qs) $.1) 
for which the hyperbolic line element transforms according to 
(ds)? = (dx)? — (dy)? — (dz)? = hi(dq,)* — hi(dq2)” — h3(dqs)* (4.2) 


where /,.2.;, are positive, real coefficients.* In the sense that the metric is diagonal, 
these transformations may be said to be orthogonal, but the absence of cross products 
like g,q; does not necessarily imply the (Euclidean) geometric orthogonality of the para- 
metric family of surfaces g; = constant with the family g; = constant. [F.g., x’ and y’ 
of (3.6) are not (Euclidean) orthogonal coordinates, but their metric does satisfy (2) 
above.| 

Introducing the transformation defined by (1) and (2) in (2.3), we have by analogy 
with Lamé’s transformation’ of Laplace’s equation 


ae nif (haha y ') _ (hale ) _ (luke | bi 
LID _— hy hehs) ( h, do Pa ( hs Pa. mS ( hs Go, " ; t.3) 


It would, of course, be possible to include 7 in the transformation, writing 
(ds)” = (dx)* — (dy)* — (dz)’ - (dr)? (4.4) 


with an obvious extension of (3) for O¢ — ¢,, . However, aside from the Lorentz trans- 
formation (2.2) we shall include 7 only in some rather simple homogeneous transforma- 
tions, where the introduction of metrical coefficients would appear rather ponderous. 


’‘P. A. Lagerstrom, Jet Propulsion Lab. Rep. 4-36; NACA, T.N. 1685 (1948). 

6‘W. D. Hayes, Thesis, Calif. Inst. Tech., Pasadena, Calif., 1947. 

*The introduction of the hyperbolic distance in the study of the wave equation is of course not new, 
having been developed in some generality by Hadamard (ref. 10, infra; also ref. 16, 430 ff.). However, the 
particular problem of separation of variables in hyperbolic space seems to have received little previous 


consideration. 
7E. T. Whittaker and G. N. Watson, Modern analysis, Macmillan Co., New York, 1948, p. 401. 








1954] THE LINEARIZED UNSTEADY SUPERSONIC FLOW 5 


5. Methods of solution. We shall concern ourselves primarily with obtaining solutions 
to (2.3), (2.4) and (2.7) by separating variables. It might be thought sufficient to seek 
solutions of (2.1), on which a considerable literature already exists, but this is not gen- 
erally the case, due principally to the difficulties associated with moving boundaries. 
Nevertheless, several interesting results may be so obtained*, and the approach has the 
advantage of physical clarity, since 7 is a direct measure of physical time, whereas r 
is not. 

An alternative attack based on the existing literature for the wave equation would 
be to rewrite (2.3) in the form 


?;; + Puy — ?:: = zz ° (5.1) 


[In view of (1), the remark in ref. 8 that there is no transformation that will fix the 
coordinate system in the wing and still yield the wave equation seems to require some 
modification.] This approach has proved quite fruitful in the steady flow case (¢,, = 9), 
due both to the physical and mathematical analogies afforded (e.g., von Karman’s 
acoustic analogy®) and, more importantly, the applicability of Hadamard’s method”’. 
(Nevertheless, in the light of subsequent developments, the most successful of the general 
methods applicable to the steady flow wing problem appear to have been those of Buse- 
mann'''’*''® and Evvard,"*’’® for which no clearly defined antecedents existed in the 
classical literature.) While Hadamard’s method is not directly applicable to the three 
dimensional wave equation"’, there exists the even more elegant method of Marcel 
Riesz’. We have not investigated the application of Riesz’s method to the unsteady 
flow problem, but some consideration has been given to this matter by P. A. Lager- 
strom’. It would appear to be of interest primarily in obtaining a solution to the direct 
wing problem (¢, specified everywhere on z = 0) and in formulating the integral equation 
for the indirect wing problem (different derivatives of ¢ specified over different parts of 
z=0Q). 

We turn now to the problem of separation of variables. By analogy with the classical 
wave equation’’’*® there must exist eleven coordinate systems in which the hyperbolic 
Helmholtz equation of (2.4) is separable. In general, these systems will differ from their 
counterparts in Euclidean space, but the cylindrical (with respect to the xz axis) systems 
remain unchanged. Since in each (cylindrical) case the separation of x yields an exponen- 
tial solution, we consider for these systems the solution of the modified equation (2.7). 


8H. Lomax, M. A. Heaslet, and F. B. Fuller, NACA, T.N. 2256 (1950). 

Th. von Karman, J. Aero Sci. 14, 373 (1947). 

10J, Hadamard, Lectures on Cauchy’s problem, Yale Univ. Press, 1923. 

14, Busemann, Luftfahrtforschung 12, 210 (1935). 

2S. Goldstein and G. N. Ward, Aero. Q. 2, 39 (1940). 

8].agerstrom, |. c. ante. 

4J. C. Evvard, NACA, T.N. 1382 (1947). 

6G. N. Ward, Q. J. Mech. and Appl. Math. 2, 136 (1949). 

6R. Courant & D. Hilbert, Methoden der mathematische Physik, J. Springer, Berlin, 1938, vol. 2, 
p. 443 

17M. Riesz, Acta Math. 81, 1-218 (1949): see also B. Baker and E. T. Copson, Huygens’ principle 
Oxford U. Press, 1950, p. 57. 

18In the unpublished lectures referred to above. 

19H. P. Robertson, Math. Ann. 98, 749 (1938). 

201,. P. Eisenhart, Ph. Rev. 45, 427 (1934); Ann. Math. 35, 284 (1934). 
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6. Cartesian coordinates. Equation (2.7) is separable in the following cylindrical 
systems: Cartesian, polar, elliptic and parabolic. The separated solutions in Cartesian 
coordinates are exponential, and the method of generalization is that of Fourier. The 
great majority of the literature on the supersonic wing problem uses Cartesian coordi- 
nates, and no further discussion is warranted here. 

7. Cylindrical polars. Let (p, ¢) be the polar coordinates defined by 


y + iz = pe'’. (7.1) 
The corresponding transformation of (2.7) yields 
p(p®,), + Poe» — (Ap)"S = 0. (7.2) 


The separated solution of (2), which differs from its classical counterpart only in the 
sign of \”, is given by 
@ = K,(Ap)e”**® (7.3) 


where K, is Macdonald’s solution to Bessel’s equation of order » and is dictated (in 
preference to alternative Bessel functions) by the boundary condition at infinity. 


A 


In the case of a rectangular wing edge (op = 0) with upper and lower surfaces g¢ = 0 
and 7, respectively, the solution 
eee ~ 
i oO - . ; a - 
- p = K..1,(Xp) Cos m+ Mt) (7.4) 
\OYy o 
is appropriate to the boundary condition 


©,” oe ie (7.5) 


(More generally, 7” can be replaced by an m-th order polynomial in y.) and may be 
used to construct a general solution for the oscillating rectangular wing. This approach 
has been used by Rott*’ for the special case m = 0, following Lamb’s treatment of 
the half plane diffraction problem.~ 
8. Elliptic cylinder coordinates. The well known transformation 
y + iz = b cosh (€ + in) (8.1) 
yields 
&,, + &,, — (Ab) *(cosh’ — — cos’ n)& = 0 (8.2) 


which separates into Mathieu equations in both — and 7. A general solution of (2), 
y; 


subject to the appropriate null condition at infinity, is given by 


@ = Gek,(é, —q)|A ce,(n, —Q@) + B,se,(n, —q)], (8.3) 
gq = 3'0'," (8.4) 
where the notation is that of McLachlan”’. 


2IN. Rott, J. Aero. Sci. 18, 775 (1951). 
2H]. Lamb, Proc. Lon. Math. Soc. (2), 4, 190 (1906); ibid 8, 422 (1910); Hydrodynamics, Dover 


Publ., New York, 1945, p. 538. 
23N. W. McLachlan, Theory and application of Mathieu functions, Oxford U. Press, 1947. 
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These coordinates have been used to obtain a general solution for an oscillating 
rectangular wing of (in principle) arbitrary aspect ratio”, following Sieger’s solution 
of the analogous diffraction problem*’. While the Laplace inversion of the resulting 
solution in terms of tabulated functions does not appear to be possible, an expansion 
in powers of \b leads to an asymptotic (in x) expansion for the potential and to integrals 
(e.g., lift and moment) that are useful for values of the “effective’”’ aspect ratio less than 
unity. [The solution for a single edge (vide supra) suffices to handle the greater than 
unity case.] 

9. Parabolic cylinder coordinates. In this case, we have 


y + tz = 3 + in)’, (9.1) 
&: + &, — VE + 1)& = 0 (9.2) 
and the resulting solution appears in the form 
® = D,{(2d)'€]D_,-1[(2d)'*n] (9.3) 
where D, is Weber’s parabolic cylinder function of order » in the notation of Whittaker”’. 
In practice the manipulation of the D, is involved, but, having introduced (é, 7), 
it may be possible to write down solutions in terms of more elementary functions. Thus, 
following Lamb (ref. 22), we find the solution 


2 (A/2)1/2(E4n) a (A/2)*472(E—y) 
is ol A + B | e* ar | +e lc + D | e* ar |, (9.4) 


0 


Evaluating A, B, C, D, by the imposition of the boundary conditions 


®.|,-92 = —W(s), y > 0 (9.5a) 

> = 0, y <0 (9.5b) 

together with the null requirement at infinity, we find for the Laplace transform of 
the potential on the upper surface (n = 0) of a rectangular wing 

|... = X'W(s) erf [(Ay)'”7] (9.6) 


in agreement with the result obtained by the method suggested in section 7 for the case 
where the boundary condition on the wing is independent of y. 
10. Hyperbolic coordinates. Let r denote the hyperbolic radius, viz. 


r= (2? —y? — 2)”, (10.1) 


Then the Euclidean transformation to spherical polar coordinates suggests the analogous 
transformation 

x =r coshé, y = rsinhé cos¢g, z = rsinhésin ¢. (10.2) 

If we restrict both r and & to be positive and real, only points within the downstream 


Mach cone [x > (y° + 2’)'”’] of the origin are included in the transformation, and the 
surfaces obtained by holding r, — or ¢ constant are circular hyperboloids (of two sheets, 


2J. W. Miles, Aero. Q. 4, 231 (1953). 
%B. Sieger, Ann. d. Physik 27, 626 (1908). 


ref. 7, sec. §16.5. 
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i?) 


although only the sheet directed along +2 is included), circular cones and planes, 
respectively. The Mach cone itself, being a characteristic surface of 0¢, has the de- 
generate specification r = 0 and = o. We remark that the tangent planes to the sur- 
faces r = constant and = constant have complementary (rather than perpendicular) 


slopes. 
Introducing the transformation (2) in (4.3) we obtain 
O¢ = (rsinh £)~*[sinh’ £(r*¢,), — sinh &(sinh é¢,), — ¢,,] (10.3) 
and, separating variables, a general solution to (2.4) is given by 
@ = 7 '°Z,41/2(0r)Bi(cosh é)e"* (10.4) 


where Z,.; is a solution to Bessel’s equation of order v + 3, and B? is a solution to the 
generalized Legendre equation 


(1 — 2)B.), + bo + DP) — wl — 2) "|B = 0. (10.5) 


Solutions of type (4) have been studied by Hayes (ref. 6) for the case of steady 
flow (» = 0), where the r dependence reduces to r’. The resulting solutions are homo- 
geneous of order », or, in the terminology of the day, “generalized conical flows’. 

An alternative attack on (2.3), after posing the dependence on é and ¢ already 
found, is to assume the solution to be homogeneous in (r/r). Thus, it is found that a 
homogeneous solution of order k is given by 


@=r ll — (2/7)? ['7°*? BE*"(7/r) Be(cosh de” (10.6) 


a result that is reminiscent of a homogeneous solution to the wave equation proposed 
by Bateman”’ and has possible application to transient loading problems. 

11. Prolate hyperboloidal coordinates. This system is derived by analogy to the 
conventional prolate spheroidal set. Modifying the transformation to the latter, we 
arrive at the new transformation 


aI 


y = ( — 1)'"(7’ — 1)'” cose 


H 94 
(EE.3) 


z= (? — 1)'"(7’ — 1)” sing. 

As in §10, the entire manifold (£, 7, ¢) is mapped in the downstream Mach cone from 
the (x, y, z) origin, but to avoid ambiguity we impose the restrictions 

l<n<é. (11.2) 
The surface of revolution § = constant, as defined by 

7 (y? +2) . 

Be. TE! am i (11.3) 
E : =~ 
is evidently a circular hyperboloid of two sheets (cf. §10) directed along the z axis. 
The same result holds for » = constant, it being necessary only to replace ~ by 7 in 
(3). Again, the respective families of surfaces are not orthogonal in the Euclidean sense. 
Introducing the transformation (1) in (4.3), we obtain 


De = & — {Ie — Dede — (Cn — Dd], — & — 1I-'(n*? — 17'S 5 6}- (11.4) 





7ref, 3, p. 384. 
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Posing the solution 


= sD glade (1.8) 
in (2.4), separation of variables yields the Lamé equation 
(@ — Dfele + [’'? — WE — "+ r)f = 0 (11.6) 


and an identical equation for g(n). 

We remark that f(£) and g(7) need not be the same solution to (6); 7.e., they may 
be Lamé functions of the first and second kinds, respectively, or vice versa. 

12. Oblate hyperboloidal coordinates. Modifying the conventional transformation 
to oblate spheroidal coordinates, we write 


x = én, 
y = # + 1)'(n’ — 1)” cos ¢, (12.1) 
gE = ( + 1)'%(n° — 1)” sing. 


In this case, symmetry exists between ¢ and in, and the surfaces § = constant and 


n = constant, v2z. 





zs _ @ +s) _ 

e-@+y7 7 
2 2) 2 

, + B m 5 - % (12.3) 


are circular hyperboloids of two and one sheets, respectively, directed along the z axis. 
The restricted range 

I<n<+@4+1)” (12.4) 
includes all points inside the downstream Mach cone, while 7 > (¢’ + 1)'” gives points 


outside of this cone. 
The transformation of the hyperbolic Laplacian yields 


De = @ + MIE + Dek — (Cn — Do), — @ + "(Cr — 1276.6} (12.5) 


and the resulting separation of (2.4) again leads to Lamé functions. 
13. Hyperboloidal coordinates. Modifying the conventional transformation to ellip- 
soidal coordinates (in the notation of Hobson”), we write 


x = h'kénf, 

y _ h-"(k? a h*) 7 =. h?)'? (9? oe hy)? —_ h’)'”?, (13.1) 

z= kk = h*) "°° —_ k?)'?(k? = "ae a k*)'” 
O<hinsksSét<. (13.2) 


The restricted range of (2) includes all points within the Mach cone. 


The coordinate surfaces are defined by 
2 2 2 
a-pip-popth (13.3) 








28, W. Hobson, Spherical and ellipsoidal harmonics, Cambr. U. Press, 1931, Ch. XI. 
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where ~ may be replaced by either 7 or ¢. The & and ¢ surfaces are elliptic hyperboloids 
of two sheets directed along the x axis (7.e., the section x = constant yields an ellipse, 
while either y or z = constant yields a hyperbola), while the 7 surfaces are elliptic hyper- 
boloids of one sheet directed along the y axis. 

The hyperbolic Laplacian is given by 


D¢ = (F -— F)'(e — 7) 'L{o} —(" -f)"€ -— wn)” 4} 
—( — Wf) '\(& — 1) 'L,{¢}, (13.4) 
Ligh = —WPy'?@ — ROE — W)'?€ — B)'deke , (13.5) 
and the resulting separation of (2.4) yields Lamé functions in each of the three variables. 
L, is identical with L; , but (¢° — k*)'”* must be replaced by (k* — 7°)'” in L, . 


14. Hyperboloido-conal coordinates. If, in (13.1) et. seg. we assume ¢ to be very 
large, so that ¢, (¢° — h”)’”’ and (¢ — k’)'’* each may be replaced by 7, we have, in analogy to 


the spheroconal coordinates of classical potential theory (see ref. 26), 


y=h'(k — Wh) ere — h’)'7 (9? — h’)'”, (14.1) 
2 (hk — Wh) 'erg — k’)'? (kh — 9)!” 
These are the hyperboloido-conal coordinates first introduced by Robinson” and applied 
by him to delta wings in both steady and unsteady flow. (Robinson uses various notations 
in the papers cited and also introduces Jacobian elliptic functions). The surfaces obtained 
by setting 7, £ or 7 constant are circular hyperboloids of two sheets directed along the 
x axis, elliptic cones directed along the x axis and elliptic cones directed along the y 
axis, respectively. In particular, the surface § = k degenerates to the triangular lamina 
bounded by y = +k '(k° — h’)'°’x and z = 0+ and lying entirely inside the Mach 
cone, thereby furnishing the desired separation of variables for the delta wing with 
subsonic leading edges. 
The hyperbolic Laplacian in these coordinates is given by 
O¢ =r [(r’¢,), — & — 1) "(Lefo} + Life}, (14.2) 
where L; and L, are given by (13.5). A solution to (2.3) that vanishes on the Mach 
cone and is regular on the z axis is given by 
@ = v,(r, 7) FOE Nn), (14.3) 
where E and F denote Lamé functions of the first and second kind in the notation of 
Hobson*’, and 
(r°y,), -nn+ly=ry,,. (14.4) 
Solutions to (4) may be obtained by comparison with the (r, 7) portions of (10.4) 
and (10.6), viz. 
(xr)e*" (14.5) 


; 


lee) see Bes 
y,(7, 7) = rl 1 — (r/r)?? ['7%*? Betz /r). (14.6) 


294. Robinson, RAE 2151, ARC 10222 (1946); J. Roy. Aero. Soc. 52, 735 (1948); Rep. 16, Cranfield 
Coll. Aero. (1948); Proc. 7th. Inter. Congr. Appl. Mech. 2, 500 (1948); cf. also Haskind and Falkovich, 
Akad. Nauk. SSSR Prikl. Mat. Mech. 11, 371 (1947) and Germain and Bader, Recherche Aero. 1949, 


3 (1949). 
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The homogeneous solutions may be applied to the solutions of the gust loading of a 
delta wing, but the practical difficulties entailed by the introduction of the Lamé func- 
tions are considerable. 

A more detailed discussion of the properties of these very useful coordinates is given 
in the papers cited*”. 

15. Parabolic coordinates. A set of coordinates closely related to the classical 
parabolic coordinates is defined by the analogous transformation [ef. (9.1)] 


c=} +7), y=éncosg, 2= ésing. (15.1) 


While the entire (£, 7, ¢) manifold is mapped inside the downstream Mach cone, 
we avoid ambiguity by imposing the restriction 


0<7<é. (15.2) 


However, the roles of £ and 7 could equally well be reversed by virtue of their symmetry 
in (1). The surfaces £ = const. and 7 = const. are non-orthogonal (in the Euclidean 
sense) families of paraboloids. 

The transformation of the hyperbolic Laplacian yields 


Lid = (¢ ™ n) ‘lé (Ede) _ n ‘(ndb,)»] = (En) “boy (15.3) 
and a general solution to (2.4) is given by 
@ = F'n Wr rpoline) Wy ro(inn e”*, (15.4) 


where W denotes a Whittaker function’. 

16. Paraboloidal coordinates. A set of coordinates that bears the same relation to 
the coordinates of §13 as the (relatively little used) paraboloidal to the ellipsoidal co- 
ordinates of classical potential theory” is given by 


c= 2VC+er +o —h — K), 
y = (k* — h’) I 2(¢? — pf)! (4? a h?)'2(¢ oo h?)'?, (16.1) 
= (k° 7 h’) : *(¢" — k’)' *(k? = n)' ci a py. 


The ranges of (£, , ¢) are specified by (13.2) for points inside the Mach cone. The co- 
ordinate surfaces are given by 


y < 


€-i ~ @-#F) 


2t— 


= &, (16.2) 
where £ may be replaced by either 7 or ¢ in (2). The surfaces £, 7 or ¢ = constant are 
respectively elliptic, hyperbolic and elliptic paraboloids directed along the zx axis. 

17. Other possibilities. There exist further coordinate systems having orthogonal 
metrics, as assumed in (4.2), but, in view of the classic investigations of Schroedinger’s 
equation (refs. 19, 20), it does not appear that separation of variables could be achieved 
in other than the eleven systems enumerated in §6-16. 

It should perhaps be pointed out that there exist coordinate systems that do not 
satisfy (4.2) but that nevertheless may be extremely useful in practice. Thus, the rotation 

%ref. 7, ch. XVI 

31J, C. Maxwell, Treatise on electricity and magnetism, Oxford Univ. Press, 1904, p. 240. 
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to ‘Mach coordinates”, (which, of course, are also Cartesian and therefore lead to solu- 
tions by separation), furnished by 


e=2“e+n, y= 2(-E+0), 252 (17.1) 

yields 
(ds)? = 2dt dn — (dz)’, (17.2) 
Od = 2¢;, — 2 - (17.3) 


These coordinates have been used to advantage by Evvard™ in attacking the unsteady 
flow problem, although it should be remarked that his end results are of questionable 
validity for time dependences other than linear.* 

In the case of steady flow (« = 0) there are many more coordinate systems in which 
(4.2) is a valid representation and O¢ = 0 can be separated. Bipolar coordinates furnish 
a cylindrical example, while toroidal coordinates (in which Laplace’s equation is sep- 
arable) can be appropriately modified. 





#J. C. Evvard, NACA, T.N. 1699 (1948). 
*In NACA, T.N. 951 (1950) it is stated that the results of T.N. 1699 are only “approximate.”’ 
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ADDITION THEOREMS FOR SPHERICAL WAVES* 


By 
BERNARD FRIEDMAN anv JOY RUSSEK 
New York University, Mathematics Research Group 


Abstract. Expansions or “addition theorems” for the spherical wave functions 
j.(kR)P™(cos @) exp (im), h\? (kR)P(cos 6) exp (im@), and h,” (kR)P(cos 6) exp (im), 
with reference to the origin O, have been obtained in terms of spherical wave functions 
with reference to the origin O’, where O’ has the coordinates (rp ,  , ¢.) with respect to O. 

1. Introduction. When a plane wave is incident on a configuration of spheres, the 
scattered field may be obtained by a series of successive approximations. The zero-th 
order approximation is just the sum of the fields scattered by each individual sphere 
when the excitation of each sphere is taken to be the original plane wave. Higher order 
approximations take into account contributions to the excitation of a particular sphere 
from the waves which the remaining spheres scatter. In particular, the first order approxi- 
mation is the field scattered by the configuration when the excitation field is taken to 
be the initial plane wave plus the zero-th order approximation to the scattered wave. 

In order to carry out this approximation it is necessary to decompose the scattered 
field of a sphere with center at O into incoming spherical waves for a second sphere 
with center at O’. Such a decomposition requires an ‘‘addition theorem” which expresses 
a spherical wave with center at O, in terms of spherical waves with center at O’. 

Such an expansion or “addition theorem” for cylindrical waves is well known [1]. 
Using this expansion, Twersky [2] was able to calculate the scattered field obtained by 
a plane wave striking a configuration of cylinders, taking into account the contributions 
to the excitations of a particular element by the radiation scattered by the remaining 
elements. The analogous case of a plane wave striking a configuration of spheres may be 
treated with the help of the addition theorems we have obtained in this paper. 

The first expansion treated in this paper is that for the standing spherical wave 
jAkR)P™(cos @)e'”® with reference to origin O. This spherical wave is represented as an 
integral of plane waves over all possible directions. These plane waves are also referred 
to the origin O. A transformation is then made to obtain the plane waves in terms of 
the origin O’. In order to evaluate the resulting integral it is necessary to have a formula 
which expresses the product of two associated Legendre functions in terms of a sum of 
associated Legendre functions. Once this is obtained by using a formula due to Infeld 
and Hull, the evaluation of the integral is possible and the final form of the expansion 
is a sum of spherical waves with reference to the origin O’. 

The second expansion treated is that for the outgoing spherical wave 
hi’ (kR)P™(cos 0)e'"® with reference to the origin O. The same procedure is used as for 
the expansion of j,(kR)P’(cos 6)e'"*. However, in this case, the integral is over plane 
waves with real and complex directions. This leads to convergence difficulties when 
the transformation from the origin O to the origin O’ is effected. These are discussed in 


*Received Feb. 11, 1953. This work was performed at Washington Square College of Arts and 
Science, New York University, and was supported in part by Contract No. AF-19(122)-463 with the 
U.S. Air Force through sponsorship of the Geophysics Research Directorate of the Air Force Cambridge 
Research Center, Air Research and Development Command. 
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Appendix III. The final form of the expansion as a sum of spherical waves with reference 
to the origin O’ is then obtained. 

From the first two expansions we then obtain an expansion for h,”’(kR)P")(cos @)¢ 
with reference to origin O, in terms of spherical waves referred to origin O’ by using the 
relationship h{'’(kKR) + h,”’ (KR) = 2j,(kR). 

We wish to express our thanks to Dr. Victor Twersky who suggested this problem. 

2. The expansion for j,(kR)P”.(cos @)e'"*. Consider a point P, which has the spherical 
coordinates (R, 6, ¢) with respect to the origin O. A plane wave coming in along the 
direction @ = a, ¢ = 8 can be expressed as follows, in terms of spherical waves with 
center at O (see [3]): 


exp (7kR cos y) = > t"(2n + 1)),(kR)P,(cos y), 


or using the addition theorem for Legendre polynomials, 


exp (2kR cos y) 


= Le » i"(2n + 1)7,(kR){n, | m |} P%(cos a)P(cos 6) exp [—im(B — ¢)], (1) 


where we have put 


(n — m)! 
(rn + m)! 
Here y is the angle between the directions (a, 8) and (@, ¢) so that cos y = cos 6 cos a + 
sin 6 sin a cos (¢ — 8). From the expansion (1), an integral representation of elementary 
spherical wave functions can be obtained. Multiplying both sides of the equation by 
P™(cos a) exp (im@) sin a, integrating over a and 8, and using the orthogonality properties 
of the Legendre functions, we obtain the well-known formula (see [3)]): 


{n, m} and P;"(cos a) = P"(cos a). 


in(KR)P™(cos 0)e'"® = (i-"/4n) | | exp (ikR cos y)P".(cos a) exp (imB) sin a da dg. (2) 
Introduce a new origin O’, where O’ has coordinates (rp , % , ¢0) with respect to O 
(Fig. 1), and let (r, 6’, 6’) be the spherical coordinates of the point P with respect to 0’. 





We shall now obtain an expansion for a standing spherical wave around O such as 
j.a(kR)P™(cos 0) exp (im@) in terms of standing spherical waves around O’ such as 
j(kr)P*(cos 6’) exp (iud’). 
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First we introduce a set of rectangular coordinates around O and a parallel set around 
O’. Let (x, y, z) be the coordinates of P relative to O, and let (2’, y’, 2’) be the coordinates 


of P relative to O’. We have 


z=R cos 8, x = Rsin 6 cos¢, y = Rsin 6sin 4, 

z’=rcos 6’, x’ = rsin 0’ cos¢, y’ = rsin 6’ sin ¢, 

g=2e’+2, z=’ +2, y=y' +m. 
Now since cos y = cos 6 cos a + sin @ sin a cos (8 — @), it follows that 


R cosy = r{cos 6’ cosa + sin 6’ cos ¢’ sin a cos 8 + sin @ sin ¢’ sin asin B] 
+ r,[cos 0 cosa + sin @ COs do Sin a cos 8B + sin 4 sin dp sin a@ sin £], 


which can be written as 
R cosy = r cosy’ + 1 COS Yo « (3) 
Here y’ is the angle between the direction (a, 8) and (6’, ¢’), while yo is the angle between 
the direction (a, 8B) and (5 , do). 
Using (3) in (2) we have 
(kR)P™(cos 6) 
| | exp (ikr cos y’) exp (ikry cos ¥o)P".(cos a) exp (im) sin a da dg. (4) 


br 


For the term exp (zkr cos y’) we use the expansion corresponding to equation (1): 
exp (ikr cosy’) = > >» i’(2v + 1){v, | w ty(kr)Pi(cos 0’)P*(cos a) exp [—iu(o’ — B)]. 


It can be shown (see Appendix I) that the above expansion is a uniformly convergent 
series and hence when it is substituted into equation (4) the order of summation and 
integration may be interchanged. Therefore we have 


7,(kR)P"(cos 6)¢ = : ys } Ee + 1){v, | uw |}7,(kr)P(cos 6’) 
eo te 


[ | fexp (tkry COS ¥o)P".(cos a)P*(cos a) exp [i(m + u)8] sin a} da as|, (5) 


In order to simplify the above equation, it is useful to have a formula expressing 
the product of two associated Legendre functions in terms of a sum of associated Le- 
gendre functions. 

Using a formula given by Infeld and Hull [4] for the integral of the product of three 
Legendre functions, we have 


P™(cos a)P}(cos a) = > ® a(| m |, | |: p, n, v)P3 "(cos a), (6) 
> 


, » — n, and where, for ‘any coefficient 


; q, n, v) we have the relationship 


(n+ v —q— I)!(2q + 1) 


Amit R GY Gag — Me +a — mig ty tnt Dl! 


where p= v+nv+tn—2,9+n — 4, >>> 
aij mi, |p 





exp [Hm +g — 9/2 + |-m| + sleil(?) in, -j — | mb, Lm +5 — ah. 


0 
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Here p = g — | m| — |u| and (s)!! = s(s — 2)(s — 4) --- 2or1, (0)!! = (—1)!! = 1, 
ng ed sé ao, 
(In Appendix II some of the a(| m |, | u |; ¢, , v) have been calculated explicitly for 
some small values of n and m.). 

Using equation (6) in equation (5) we have 


a 6)e""® 
>> 8 p> {i’(Qv + 1){v, | uw |}a(| m |, | w |; p, n, »)j,(kr)P%(cos 6’)e***' I") (8) 


where J?" is an integral of the form 


a2e ° 


Io? = | | exp (ikry cos yo)P2*"(cos a) exp [t(u + m)8] sin a da dp. 


It is possible to evaluate these integrals directly [5] (compare equation 4), and we have 
Ip" = 4nij,(kr.)P>* "(cos 6) exp [1(m + p)do]. (9) 


Substituting equation (9) into (8) we have for the expansion of a spherical wave 
with respect to one origin, in terms of the spherical waves with respect to another origin, 
the formula 


jAKR)P(cos 6)e"* = > > D le" + Diy, | [hal] m |, | a lem, m,») 
Jolkro)j (kr)P5*"(cos 6)P%(cos 6’) exp [i(m + u)do] exp (—tud’)}. (10) 


3. The expansion of h{’ (kR)P%(cos 6)e'"* and hi? (kR)P™(cos 6)e'"*. In order to 
obtain an addition theorem for the outgoing spherical wave h,'’(kR)P%(cos 6) exp (im¢) 
we must first get the integral representation for this function. We have the formula [6] 


ep e/2—-1@ 


2r 
exp (tkr,)/ikr, = es exp (tkr, cos y,) sin a da dB, (11) 
2r /0 /0 


where 7; , 9; , ¢; are the coordinates of a point P, with respect to the origin O, . Now, 
let the point P, have spherical coordinates (R, 6, ¢) with respect to the origin O, and 
let O, have coordinates r°, ¢°, 6° with respect to O. From equation (3) we see that 

R cosy =r cosy’ +7, cosy, 


or 
exp (ikr, cos y,) = exp (—ikr’ cos y’) exp (ikR cos 7) (12) 


where as before cos y = cos a cos @ + sin a sin 6 cos (8 — ¢). 
Now from equation (1) we have 


exp (—ikr° cos ’) 


= >> Dd i"(2n + 1){n, | m |}j.(kr)P'2(cos 6°)P™(cos a) exp [im(6 — ¢’)]. (13) 


n=0 m=—n 


Therefore, using the formula [7] 
exp (ikr,) /ikr, 
—— > (2n + 1){n, | m |}j,(kr’)h,”(kR)P2(cos 6°)P™(cos 6) exp [im(@ — ¢”)] (14) 


n=1 m=—nfn 
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and equation (12) and (13) in equation (11) and comparing the coefficients of each | 
j.(kr°)P™(cos 6°) we find the integral representation to be [8] 


hy” (kR)P7(cos 6) exp (imd) 
n p~2e a e/2—-i 


-i | | 
2m . 


Now let the point P have spherical coordinates (r’, 6’, ¢’) with reference to another 
origin O’, where O’ has coordinates (rp , 9% , 0) with respect to O. (See Fig. 1). From 
equations (1) and (3) we see that 


‘ exp (ikR cos y)P(cos a) exp (im) sin a da dg. (15) 


R cosy = 7r cosy’ + 7% COS Yo 
and that 


exp (ikr cos y’) 


= 5D + df, | w [19.(br)P4(cos 6)P*(cos a) exp [—inle’ — A]. (16) 


When these two relations are substituted into equation (15), it can be shown (see Appen- 
dix III) that it is possible to exchange the order of summation and integration provided 
that r < r, . Therefore we have 

h,’ (kR)P™(cos @) exp (im¢) = _ > +B {i Qv + 1){v, | uw |}7,(kr)P(cos 6’) 


vy=0 p=-? 
r/2-i@ 


Qr 
emp (—in) [ [ exp (ikro cos 6,)P#(cos a)P"(cos a) 


-exp [7(u + m)8] sin a da as} for P<. (17) 


Now we have the formula (see equation (6)) 


P*(cos a)P™(cos a) = >> a(| m |, | u |; p, n, »P2*"(cos a), 
i 


where p = vy +n,» +n — 2, +--+,» — mand the a(| m || « |; p, n, v) are given by 
equation (7), substituting the above expression into equation (17) we have 


hi’ (kR)P™(cos 6) exp (im) = — + XS Kr + If», | u |}al| m |, | w |; p, 2, ») 


* j.(kr)P;(cos 6’) exp (—ing’)K,"" (18) 


where the K>’” are integrals of the form 


a2 a e/2-1i@ 


Kem = | | exp (ikrp cos 6,)P**"(cos a) exp [i(u + m)6] sin a da dp. 


“0 


These integrals can be evaluated directly. Using equation (15) we find that 
Ké*" = 2wi”hs” (kr,)P2*™" exp [i(m + p)do]- 


Therefore, equation (18) becomes 


hi? (kR)P™(cos @) exp (img) = ¥" > Dd {a"*?-"(Qv + 1) fv, | w |fa(| m|, | w |; p, 2, ») 


7=0 w=? P 


-j,(kr)hs” (kro)PS*"(cos 0)P%(cos 6’) exp [i(m + u)bo — tud’]}, o EM. 
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When r, < r, the resulting expansion would be the above equation with r and r, inter- 
changed. Therefore we have 
h?(kR)P%(cos 6) exp (im) = SY DY 2?" + Df, | uw [Hal ml, lw lip,n») 
y=0 wp=-? Dp 
- 7,(kr-)hs” (kr,)P2*"(cos 6)P3(cos 6’) exp [i(m + u)dbo] exp (—iud’)}, (19) 
where r- = minimum (r, 7) r; = maximum (r, ro). This is the expansion for a spherical 
wave hh,’ (kR)P™(cos 0) exp (im¢) about one origin, in terms of spherical waves around 
another origin. 
To find the expansion for h{” (kR)P™(cos 6) exp (im@) we use the relationship 
AP (RR) + hi? (kR) = 2j,(kR) (20) 


and the expansions for j,(kR)P™(cos 6) exp (im@) and h,'’(kR)P(cos 6) exp (im@). 
This gives us finally 
hy? (kR)P™(cos 0) exp (im¢d) = 2. > Zz, fy, | w |}(Qv + la(| m|, | uw |; p, 2, ») 

y=0 p=-P Dp 


- 7(kre)he (kr,)P2*"(cos 0.)P(cos 6’) exp [7(m + u)do] exp (—iud’). (21) 


APPENDIX I. 


To show that the expansion 


exp (ikr cosy’) = >) >> 1’(2v + 1){», | u |}j,(kr)P2(cos 6’)P*(cos a) exp [—in(¢’ — 8)] 
; (22) 


is a uniformly convergent series, we first write it in the form 
exp (ikr cos y’) = >. i’(2v + 1)j,(kr)P,(cos 7’). (23) 
v=0 


Now we know that for 0 < y’ < 22, | P, (cosy’) | < 1, and also| exp [—zu(¢’ — 8)] | = 1. 
For large values of », | j,(kr) | behaves like (kr/2)’/v!. Hence we have 


> | i’(2v + 1)j,(kr)P,(cos y’) | < Do (2v + 1)(kr/2)’ /r! (24) 


r=0 vy=0 


However, >_5-, (2v + 1)(kr/2)’/v! is a convergent series, and therefore it follows that 
> i’(2v + 1)j,(kn)P,(cos 7’) (25) 


is a uniformly convergent series. 


ApPEenpDix II. 

For general values of v, n, | u | and | m |, a formula has been obtained (see equation 7) 
for the coefficients a(| m |, | u |; p, n, v). Using this formula, we have calculated the 
coefficients for some special values of n and m. It is sufficient to evaluate the 
a(| m |, | u |; p, , v) for positive values of m only, since from these the coefficients for 
negative values of m are readily obtainable and depend upon the definition of P, (zx) 


in terms of P(x). 
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We see that for a(| m |, | u |; p, n, vy) = O unless p = » + lory — 1: 
























































n=1 | a(| m |, |#|;» +1, 1,») a(| m |, |#|;» — 1, 1,») 
m=0 | @—lel+D/Q%+)) + [a D/ +1) 
m= 1 | 1/(2v + 1) —1/(2» + 1). 
For n = 2, the a(| m |, | uw |; p, 2, v) = O unless p = v + 2, » — 2: 
n=2| a(|m|,|u\|;» + 2, 2,») a(| m |, | |;¥,2,»)} a(| m |, | uw |;v — 2, 2, ») 
mn = OO eal HL + De = |v] +) OG +y—3 [ne /) Be+ |e Dot |o1|—-V 
2 (2v + 3)(2» + 1) (2v + 3)(2» — 1) |2 (2v + 1)(2» — 1) 
a @=|el+D 3 v+{u) _3__@+|n) 
” (2v + 8)(2» + 1) 2 (2v + 3)(2» — 1) (2v + 1)(2» — 1) 
_ 4 3 ~6 3 
om Qv + 3)(2» + 1) (Qv + 3)(2» — 1) (Qy + 1)(2» — 1)° 








When n = 3, a(| m |, | u |; p, 3, v) = Ounless p = » + 3,9 +1,» — 1,» — 3. 
For m = 0 


- _5¢—|4|+30— lu |+20— |e#(|4+D 
lh atliaeiaate ( + DQ + 3) +1) 


— _ 8e¢-—|el|+ DC’? +2 — 5/4 /’) 
a(|m |, luli» + 1,3, %) = 9° To, 4 By + DQ — 1) 


3+ |» |e? — 5 |u|? — 1) 
(2y + 3)(2y + 1I(2v — 3) 


v+iePo+ |u| — DO+/e| — 2) 
(2v + 1)(2» — 1)(2v — 3) ‘ 











Il 


a(| m|,|u\|;» — 1, 3, ») 


bo 


feb | 





a(| m |, |u |;” — 3, 3,») = 


bo] 


AppEenpIx III. 


In order to prove the validity of exchanging the order of summation and integration 
in equation (17), we first write it in the form 


“—n 2 a«x/2-iC 


hy’ (kR)P'x(cos 0) exp (im¢) = = lim | | 


‘ 
2 com 0 “0 


{[ > i’(2v + 1)7,(kR)P,(cos ) | exp (ikro COS Yo)P":(cos a) exp (im§) sin a da as}. (26) 


However, | j,(kr) | < (kr/2)”/v! and | P,(cos y’) | < | cosy’ + (cos’ y’ — 1)” |’; 
hence 
ee) @ 9 , 2 _ 1/2 |2 
| #2 + DIER)P,(cosy) | < 5 + D Le cosy’ + (cos y” — 1) | 


r=0 y=0 y!2’ 
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where the right hand side is a convergent series, and therefore 
> i’(2v + 1)j,(kr)P,(cos y’) 
is a uniformly convergent series for any finite values of r and y’. 
Now, in the finite region containing the segment of the contour of integration in- 
cluded between 0 and 2/2 — iC, exp (ikry cos yo)P™ (cos a) sin a exp (im) is bounded, 
and hence the series 


exp (zkry cos yo)P".(cos a) sin a exp (ims) =: t’(2v + 1)),(kr)P,(cos y’) 
is uniformly convergent in this region. Therefore we can write equation (26) as 


h,”(kR)P™(cos 6) exp (im¢) 


=. — lim LS i"(2v + 1)j,(kr) | P,(cos y’)P",(COs a) 


v=0 “0 *0 


) 
-exp (ikry COS Yo) exp (tmB) sin a da dB?. (27) 
) 
Equation (27) can then be put in the form 
hi? (kKR)P™(cos 6) exp (imo) = = >> 1"(2v )7 (kr) | | P,(cos y’)P".( cos a) 


‘exp (zkro COS Yo) exp (tm) sin a da dB — S- lim 1 >. i’(2v + 1)j,(kr) 
“7 CHa v= 


| | P,(cos y')P’.(cos a) exp (ikry cos ¥o) exp (im) sin a da dg. (28) 
Jo ) 


©? #/2-iC 


If we can now show that 


K==— — lim 4 + i’(2v + 1)j,(kr) 


a2e ap e/2-1i@2 


. ° | P,(cos y’)P”.(cos a) exp (ikrp COs yo) exp (im) sin a da as} = 0, (29) 


“0 / #/2-ic 


then the proof will be completed. 
To prove that K — 0 as C — ©, we first substitute into equation (29) the expansion 


P,(cos y’) = > a fv, | w |}-P%(cos 6’)P%(cos a) exp [—ipn(¢’ — B)]. 
(Since this is a finite sum, we may interchange the order of summation and integration.) 
This gives us 


K =— — lim {> a *(2v + 1){v, | uw |}7,(kr)P*(cos 0’) exp (—iyud’) 
e2r np F/2-i@ \ 
. ° | P*(cos a)P™(cos a) exp (ikry Cos yo) exp [i(m + y)B] sin a da dB. (30) 


“0 ¥ #/2-iC 
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Now, in the (a, 8) plane we perform a rotation such that 
cos a’ = COSY) = COSa COs 6 + sin asin 0 cos (B — do) 
B’ =B-—d. 


Then cos a = cos a’ cos 6 + sin a’ sin @ cos (8’ + ¢) and the limits of integration 
remain the same. Hence K becomes 


K : _ 


- lim {s > i"(Qv + If», | w |}j,(kr) exp (—ind’)P4(cos 6’) 


| | [exp (zkrp cos a’)P*(cos a’ cos 6 + sin a’ sin 6 cos (8’ + ¢o)) 


/ r/2-i€ 


-P™(cos a’ cos 6) + sin a’ sin 6, cos (B’ + @o)) 


-exp [7(m + u)(B’ — d)] sin a’ da’ ast. (31) 


Now for large values of | z |, we have the following asymptotic form for P$(z), [10] 


pry) w 2T@ + 1/2) 
”  (p— | q|)'rd/2)’ 
If we let a’ = 2/2 — iv, where C < WV < ©», then for very large ¥, cos a’ behaves 
like ie*/2 and sin a’ behaves like e*/2. Hence, using the above asymptotic form we 
find that 








| o'’TP / ¥/92 v | 
P*(cos a’ cos 6) + sin a’ sin 6, cos (B’ + ¢o)) | ~ | fs v 7 id ia 
= | gs 177 /é 


2"T(n + 1/2)(e*/2)” | 
| (n — | m |)!T(1/2) 








P™(cos a’ cos 6 + sin a’ sin 6) cos (8’ + 6)) | ~ 


Therefore, on substituting the above expressions into equation (29), and then taking 
the absolute value of K and integrating over 6’, we have 
t"T(n + 1/2) 


K |< QnT(1/2)T 2)(n — | m \)! | 


( . » “¥/« . ‘ | 
stim} > |! = eT: 1/2) | (er)P*(cos 6’) exp (— ind’) | 3 


, | (e*)’*" exp (—e*kr,/2)e* aw}, (32) 
The integral in equation (30) is now a portion of the real axis. Letting ¢ = e* the 
integral becomes 


[ t** exp (—e"kro/2)e dv = [| '** exp (—thro/2) dt (33) 
“D 


“D 
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where D = e° and therefore as C > ~, D — o. Integrating equation (33) by parts, 
it becomes 


/ t’*" exp (—tkr,/2) dt 
D 





t’*"~* exp (—tkr,/2) as. (34) 








ae a a 
om { 2 sn 6 tkro/2) = 20+ n) [oe 
Therefore 
[ (%)""* exp (—e"kra/2)e" ay < 22" op {= Dire/2) (35) 
and using this inequality with equation (32), the result is 
| K | ie + 1/3 lim D” exp (—kr,D/2) 


|S Fara 2)(n — | m ie 
See A liu(hr) | | Pe(cos 6”) | D’. (86) 











| au |)! 
However, 
DH ia v + ae 
| Pr(cos 6) | S x ’*(2 sin 6 yw T(v + 3/2) 
and 
e: (kr /2)” 
| jkr) | a SEED 


for large values of v. Putting these values for | P?(cos 6’) | and | j,(kr) | into equation 
(34) and writing (v + 1/2)I'(v + 1/2) as ['(v + 3/2), we have 


IK|< za Pn + 1/2) —; lim D" exp (—kr.D/2) 5 ene 2 ; (37) 


= x kr,(2 sin 6’)'7(n — | m |)! pe | 





However 


o & 9 a 
a = exp (krD/2), 





y=0 


and therefore equation (37) becomes 





|K |< A lim D" exp | Dir — r)k/2 | (38) 
D2 
where A is the constant 
(n 7 ae. a 
ews sin 6’)'“*(n — | m |)!" 


If we now let D — =~ we see that | K | — 0 provided that r < r, . Therefore equation 
(28) becomes 





hy” (kR)P'(cos 6) exp (im@) = )j(kr) 


27F 2/2-i@2 
. [ | P(cos y’)P(cos a) exp (ikro cos 70) exp (imB) sin a da d8 whenr < ry (39) 


70 “0 


and when P,(cos y’) is expanded in terms of P?(cos a)P?(cos 6’) we have equation (16). 
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TWO DIMENSIONAL SOURCE FLOW OF A VISCOUS FLUID* 


BY 
H. C. LEVEY 


Aeronautical Research Laboratories, Melbourne 


Summary. The steady two-dimensional source-type flow of a viscous heat-conduct- 
ing perfect gas is investigated. The solutions of physical significance all contain shocks, 
and bounds are given for the shock-thickness in terms of the shock-strength and the 

teynolds number of the flow. It is found that the entropy rises to a maximum within 
the shock, and this maximum does not disappear even when the viscosity tends to zero. 

Introduction. The main part of this investigation is concerned with two-dimensional 
source flow similar to the flow in a divergent channel with straight walls, for instance, 
in a nozzle or a diffuser, when boundary layers are neglected. If the fluid is assumed 
inviscid, no fundamental length exists, but if it is viscous a Reynolds number R char- 
acterizing the flow is provided by the ratio of the mass flow (per unit length normal 
to the plane of the flow) to the viscosity (cf. eq. (1.14)). In ordinary supersonic wind- 
tunnel nozzles this Reynolds number is of the order of 10’, but in low-density, hyper- 
sonic wind tunnel nozzles (for which indeed the conical shape is being increasingly 
favoured) it can be of very much smaller order, and deviations of the flow from that 
obtained in the limit R — © may be of some interest. The problem has also been con- 
sidered by Sakurai [1], who derived an equation similar to (1.41), and sketched the 
solution curves for R =~ 20 (see 1.4). 

The problem is, moreover, of theoretical interest since the corresponding problem 
for an inviscid fluid is one of the few for which an exact solution containing a limit 
line has been found [2]. This limit line is the sonic circle, and its exterior is doubly covered 
by the velocity field (1.2). The subsonic branch of the solution has a stagnation point 
at infinity and we aim to find the corresponding solutions for a viscous fluid. 

The energy equation is integrated once to give two first order simultaneous differ- 
ential equations for V and @ as functions of w, where V is essentially the velocity gradient, 
and w and @ are respectively the speed and temperature in dimensionless notation 
(1.1 and 1.3). All the solutions which have a stagnation point at infinity are shown 
to have formal asymptotic expansions for V and @, in this neighbourhood, which agree 
to the first order with the inviscid solution when R — @ (1.3). 

The simultaneous differential equations for V and @ are of the singular perturbation 
type. The highest derivatives are multiplied by a small parameter, namely R™', and 
for a first approximation to equations of this type, the small parameter is taken to be 
zero. It is fairly obvious that as long as the highest derivatives remain ‘small’, there 
are solutions of the full equations which differ little from the solutions of the lower- 
order approximate equations. But it cannot be expected that the full solutions will 
converge uniformly to the approximate solutions over the whole flow field, as the param- 
eter tends to zero. The boundary layer is a case in point, and in general, one may expect 
regions where the limiting solution is quite different from the inviscid approximation. 

For a particular value of the Prandtl number o, the energy equation can be inte- 
grated a second time, and the simultaneous equations are then reduced to one first 


*Received April 4, 1952; revised manuscript received March 4, 1953. 








26 H. C. LEVEY [Vol. XII, No. I 


order, non-linear differential equation for V in terms of w, which is still of the perturbed 
type (1.4 and eq. 6.5). This value of o is 3/4 + 0(R™), and as o = 0.72 for air, the 
results obtained below should be significant. 

In sec. 2, 3 and 4, the solution curves of this differential equation are discussed in 
the w-V plane by semi-geometrical methods for fixed, large R. Simpler curves are 
given which provide bounds enclosing the solution curves. In the limit as R — ~, the 
singularity found is not a limit line, but an ideal shock joining parts of inviscid solutions. 
A limit line does not occur even for the exceptional solutions, physically unrealisable, 
which do not contain a shock. 

A full discussion is given for the case of constant viscosity. For the case of a viscosity 
proportional to the absolute temperature, the method of solution is outlined, and the 
results are summarised in sec. 6. In practice the variation of viscosity with temperature 
lies between these two extremes. There is no qualitative difference between the two cases 
as regards the shock formation. However, in the case of variable viscosity, the absolute 
temperature is automatically positive throughout the flow if it is positive at any one point 
of it, whereas in the case of constant viscosity there are solutions violating this require- 
ment which have to be discarded on physical grounds. 

1. The fundamental equations. 

1.1. In this section, the equations for the steady flow of a perfect gas are first con- 
sidered in general, in order to derive a first integral of the energy equation. The equations 
are then specialised to the case of purely radial, two-dimensional flow. 

Let 27; , X2 , 2; be a right-handed system of Cartesian coordinates, let v; be the com- 
ponent of velocity in the direction of z; increasing, and let p, p, T, u, A, R, C, , C, and 
a, denote respectively the pressure, density, absolute temperature, viscosity, heat- 
conductivity, gas constant, specific heat at constant volume, specific heat at constant 
pressure, and local speed of sound ((dp/dp)s). The equations for the conservation of 
mass, momentum and energy are, with the usual summation convention, 








*) 
— (w;) = 0, (1.1) 
Ox; 
ov Fi) ( 2 2s) rs) 
y= - = ep ——16;; 2 —— (pei), 1.2 
ie Ox; Ox; P+ 3 - OX, ai Ox; (ues) (1.2) 
and ; : 
a, Ov, j 2 av; dv; \ Fi) ( af) 
nan CF a me At =e oT + A ah Lf 
al Ox; ( 1) + P OX, aa ate 3 Ox; oH + Ox; Ox; ( 3) 
where 
1 (2 2) 
3 = 7 —}. 1.4) 
Ci 2 \dz; + Ox; (1.4, 


and 4,; is the Kronecker delta. Multiplication of (1.2) by v; and addition to (1.3) gives 


2 | wiC.T + & pod 2 yo, — awe — 022} = 
ax; {mn C.0 + 2 pu v; + PV; + 3 Mv; ax, 2Quv €;; X ax, = 0, 


and so 





i 2, i. ve. seek f 
pu.C,T + 2 pu v; + PN; + 3 MU; Oz; Qyv Ci; A Ox; = sik Ox; ’ (1.5) 


7 


where the A; are unknown functions of the 7; . 
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The flows investigated in the following are those which are independent of z; , and 
in the (x, , x2)-plane, depend only on the distance, r, from the origin. The only non- 
vanishing velocity component is that in the direction of r increasing, denoted by uw. 
Referred to polar coordinates, (r, ¢) in the (x, , 22) plane, (1.5) then yields 





pe 1 2) 2 ( “) du _ of _ 190A; 
mcr +2 +hu ae a eo ee ir 
and 
_ As 
0= = (1.7) 


Hence A; is a function of ¢ only, and 0A;/dg must be a constant. 
From equation (1.1) 
pur = «x (a constant), (1.8) 
and (1.2) gives 


, ap _ 44 / au) as de _ te _ 2.4 (em) (1.9) 
or ' dr 38dr\" ar r dr r S38ar\r/° ; 
With the help of the equation of state of the gas, 
p = RoT, (1.10) 
(1.6) may be written 
12) 4 2afe M4) — auur  — eo 
pu C,P 5 oo 3 Mu\T +u Quur ar Ar Y “a v. (1.11) 


Equations (1.8) to (1.11) govern the steady, two-dimensional, purely radial flow of a 
perfect gas. 
These equations are put into non-dimensional form by the substitutions 
ot a" ue (r+4)": —/agaa _2 - 
w= ( 9 —" 2 » 25° 6= T,’ & = logr. (1.12) 


The point at infinity will be taken to be a stagnation point, and 7, and a» denote the 
corresponding temperature and speed of sound, respectively; y = C,/C, . The Prandtl 
number ¢ is defined by 





ie wt (1.13) 
and the only remaining dimensionless parameter is 
38ky +1 
Ra +, : 
a dy (1.14) 


which may be regarded as the Reynolds number of the flow. 
Elimination of p and p from equations (1.8) to (1.11), and use of (1.12) to (1.14) 
leads to the two equations 


(8 + 1)w*w’ + wl’ — wd — bw’ = & w'(w!" — w) + 2w*(2w’ — w) - (4), (1.15) 


—— is ” 3 "en 88 
ia ot + Ra + me RA + Bo’ ~ RA +8) 





ww’ = C (aconstant), (1.16) 
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where 


B=y7-)D)/ot+D (1.17) 


and a dash denotes differentiation with respect to &. 

1.2. The “‘inviscid solution.” The solution for source flow in an inviscid perfect gas 
is obtained by putting R = © in the above equations, without any enquiry as to the 
validity of such a step. It will be seen later that some of the terms which have a factor 
of order R~’, are themselves large of order R, so that even when R becomes infinite a 
finite contribution remains. If this latter possibility is ignored, the equations reduce to 


w{(1 + B)ww’ + 0’ — 6} — Aw’ = 0 (1.18) 
and 
6+ Bw = 1, (1.19) 
since 6 = 1 at the stagnation point at infinity by definition. 
When @ is eliminated from (1.18) and (1.19) the equation 
(1 — Bw’) 
y= —wi— 1.20) 
u u (I — w’) ( 
is obtained. Its solution is 
: ra 1/2 . 
page MoE .g ~ ayer. (1.21) 
Ap PoW 
where p, is the stagnation density at r = @. 


The graph of w(r) is given in Fig. 1. The ‘sonic speed’ (which is the fluid speed at 


Ww 


Bt 











| a 
N Tr 
Fic. 1. w versus r for inviscid gas. 
which the fluid speed equals the local speed of sound) corresponds to w = 1, and the 


maximum speed attainable, at which the temperature falls to absolute zero, corresponds 
tow = 6”. If r, is the value of r for which w = 1, (1.20) gives no solution for r < 1, , 
but for r > r, there are two branches of the w — r curve, one tending to zero, and the 
other tending to w = 6°” as r tends to infinity. Thus for any r greater than r, there 
are two possible values of w, one representing a supersonic speed, and the other a sub- 
sonic speed. The streamlines are radial and, in fact, are cusped at the sonic circle, which 
is a limit line (of a rather special type). 
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It is intended now to use this solution as a guide to a study of the solutions of the 
equations (1.15) and (1.16), inasmuch as it is expected that sourceflow, even in a viscous, 
heat-conducting fluid, should lead to a stagnation point at r = ©, and that viscous 
effects should become comparatively unimportant there. But a limit line cannot exist 
in the real fluid, and so, if solutions of (1.15) and (1.16) are sought which behave like 
the solution of (1.20) for large r and small w, their continuation backward in r should 
show what phenomena are to be expected in a real fluid when inviscid theory predicts 
a limit line. To preserve the correspondence between the solutions for large r, the con- 
stant C in (1.16) is taken to be unity. 

1.3. The solution at large distances from the source, for large, but finite, values of R. 
In this section, » is taken to be a constant, which leads to a considerable simplification 
at this stage. When the complete solution is discussed later, however, the case of u 
varying with temperature will also be considered. It should be noticed here that the 
ratio of the specific heats, y, is assumed constant throughout, and the Prandtl number, 
o, is assumed constant when the viscosity is constant. 

It is possible to eliminate any explicit dependence on ~ from the equations (1.15) 
and (1.16) by the substitution 


,—_ 9 Ww _ _o, dw 
J 2 2r Tr? (1.22) 
which permits us to reduce these equations to two simultaneous first order equations 
in V, 6 and w, namely 
b cet 1 ied See Se . ae 
pv J _* — 51 + B)w'I — 9 Vw, wat 5 Ot + , (1.23) 
and 
g@—1+ a +s 2 rw + — a wV =0. (1.24) 
Ve RG + B)) 2Ro(1 + 8) 


The variable V is closely related to the fluid acceleration, and satisfies a certain non- 
linear differential equation which is obtained by eliminating @ from (1.23) and (1.24). 
The equations (1.23) and (1.24) respectively may be put into the forms 


‘Pe ee So 
ys (V — 2w) ow? V (I 2w) 
__R§ a ae ee 
2 yi +8 aoe t t+ be 1) 
(1.25) 
1+24\w Rj, 3 £\dé 
+ rs ‘tae a} y ~ te {1 a+ was ie 
,, 4 
i s( » 8, dw’ w) 
and 
Pe =o gy i us gee! Bl mags 
ap (F —1)+ 3) (1 + B)(0 l) = a (1 + Ay + RG . he 3 
(1.26) 


= g(V, w). 
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If z(w), ¢(w) are defined by 


" 
zw) = — | mV dw, (1.27) 
2Re 

(uw) = [ 2R2 1 + 6) du, (1.28) 

then equations (1.25) and (1.26) have the formal solutions 
V — 2w = otf etf-aw + Cie* (1.29) 

and 

6—1l=e° [ eg dw + Cie". (1.30) 


These integral equations may be solved iteratively in the neighbourhood of w = 
V = 0, starting with the first approximation 


V = 2w (1.31) 


suggested by the inviscid solution. If (@ — 1) is to be finite near the stagnation point, 
C, must be zero, and the first approximation to (1.30) is 


12 + RU + 8) — 
6 + Ro(l + 8) - 


after substitution of this expression in f(v, 6, w), (1.29) gives the second approximation 


(1 — 6’)Ro + 6(1 + B) — 248c\| , = os 
= 2f( (1+ BRo + 6 20 + Cie’. (1.33) 


This process can be repeated as many times as desired to obtain higher approxi- 
mations, and gives the asymptotic solution for V and @ near w = V = 0. The coefficients 
of the powers of w in the series so obtained agree with those found for the ‘inviscid’ 
solution, except for terms of order R~*. However, the former series are divergent, as 
the term of order R™ in the coefficient of w* is unbounded as n — . Also there is an 
infinite number of solutions having the same sort of behaviour near w = V = O, due 
to the presence of the term C,e™* in (1.33), where C, is arbitrary. This is of course due 
to the fact that the differential equations leading to these expansions are of higher 
order than the ‘inviscid’ differential equations. The term C,e* is asymptotically of 


the form 





6—1= —£e 





V — 2w 


Ci exp {- 5 Rw - ; (1 + A)R log u}, 


and so is very small if w is small enough, whatever the value of C, . 

Thus there is an infinite number of solutions for the flow of the viscous heat-con- 
ducting gas, which, in the neighbourhood of the stagnation point at infinity, are asymp- 
totically equal to the solution for the inviscid gas when the viscosity and the heat- 
conductivity tend to zero. It is very difficult to investigate the solutions of (1.23) and 
(1.24) in other parts of the flow field, but it is found that there is a particular value of 
o which enables the problem to be reduced to that of the solution of a first order non- 
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linear differential equation. This equation is still of the singular perturbation type, 
but its solutions can be discussed, after some trouble, and it is with this problem that 
we shall be concerned hereafter. 

The equation (1.24) may be written 


tes an ee Se a 
nae + Al T RG ‘hu G+ prea? --RA+DE r= 0, (1.34) 


and this equation, as it stands, can be integrated when 





o see WEE 
¢= 3 4 RA +8)" (1.35) 


For, when £ is defined by 


E=6éo-1+ ad + (1.36) 


4 2 

RO + aS ve 
the equation takes the form 

’ 4 _ dE _ - 

"~-mantia’™ (1.37) 


which has the general solution 


E = A exp [2 (1+ @)R+ |} (1.38) 


where A is an arbitrary constant. Now, E must tend to zero at the stagnation point 
at infinity, and so A = 0 and £ = O. In fact, E£ is effectively the difference between 
the total energy, per unit mass, of the fluid and the value of this total energy at the 
stagnation point, and (1.38) shows that K — + © at this stagnation point unless F = 0. 
Therefore, when o has the value given by (1.35), 


This permits us to reduce the problem to a first order differential equation for V as a 
function of w. 

The value chosen for ¢ is actually not far from the truth, as experimental values for 
air are in the neighbourhood of 0.72, and when R is large, the value chosen will be close 
to 0.75. From the form of (1.15) and (1.16), (with constant FR) it follows from [5] that 
the solutions are continuous in ¢, and in fact have continuous derivatives of any order 
with respect tog inw > 0, —©° <o < ©. In the neighbourhood of w = 0 the variation 
with o can be seen from (1.33). Thus the only effect of choosing the above particular 
value of o is to simplify the equations rather than materially affect their solutions, and 
we expect that the solutions investigated hereafter will be representative of—and, in 
fact, very close to—the actual flows which would occur in air under the same boundary 
conditions. In [1], Sakurai fits values of o) to experimental values for various gases by 
choosing appropriate values of R, (negative for air), and sketches the solution curves 
for oo = 0.88 (R = 20), Meyer’s theoretical value. 

With the help of equation (1.39), 6 and its derivatives are eliminated from equation 


(1.23) to give 


weer bode 2 ane } wil ( 41 + 26) ‘| 
R 5M (i+ RQ + 8) B+pige/@s (1-40) 


du 
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For convenience in the subsequent algebra, terms which are genuinely of order 
R~ compared with terms of order 1, are now neglected on the right-hand side of (1.40). 
That is, in each of the curly brackets, the coefficients of w* are replaced by 1 and B 
respectively. Thus, the equation to be investigated is 

ae uo = ; Vil — w) — w(l — Bw). (1.41) 
The methods hereafter applied to (1.41) are quite applicable to (1.40), but the algebra 
becomes much more involved, merely because of the unwieldy forms of the two coefficients 
of w’. The difference between the two equations is very small for large R, and physically 
the difference is that in equation (1.40) the sonic speed and maximum speed vary slightly 
with FR, while in equation (1.41) they are fixed at their inviscid values. The variation 
in w between the two is only of order R~', while in the subsequent discussion the most 
significant variations in w are of order R™’”’. 
2. Properties of the solution curves. 

Although equation (1.41) is of a simple appearance, it cannot be integrated to find 
closed solutions, so that the discussion requires indirect methods. In the first place, 
in order to investigate the form of the solutions we use the standard techniques of 
curve sketching to discuss the behaviour of the solution curves in the w-V plane, and 
to this end the behaviour of the curves of zero slope and zero curvature is established 
in 2.1, while 2.2 gives some of the more elementary properties of the solution curves. 

2.1. Let C; be the curve on which dV /dw is zero in equation (1.41), that is, the curve 
given by 

V = f(w) = 2w ae = : (2.1) 
: l— w 
which is also the inviscid solution curve (1.20). 

Let C, be the curve on which d’V/dw” is zero in equation (1.41), that is, the curve 

given by 


V* — w(1 + Bw’) V’* — 4 R11 — w’)(1 — Bu’) V + Ru(l — Bu’)? = 0. (2.2) 
It has the following properties. 


(i) When w = 0 (stagnation point), V = 0 or +(R/2)'”’. 
(ii) When w = 1 (sonic speed), 


V = -—(1 — 6)*°R'?* {1+ 0(R)}. (2.3) 


(iii) When w = 6°” (maximum speed), V? = 0 and V = 26°”. At V* = 0, there 
is a double point, the curve having there the slopes —3(1 — 6)R{1 + O(R™)} and 
4B6/(1 — 6){1 + O(R™")}. 

(iv) For w > R’” there are three branches of the curve, given asymptotically by 


V ~ 26u, 


V~- 


bol 


Rw, 


V ~ Bu". 
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(v) The curve C, has infinite slope at (w,, , V,,) where 
WwW, ™ 1 — 5h (1 — a, 


V. tas ats in FT. (2.4) 
(vi) Provided that 
[1 —w? |'R" = o(1) (2.5) 


there is a branch of C, given by 

V = g,(w) = f(w){1 + O(R"*(1 — w’)~*)} (2.6) 
which lies close to C; in the range 0 < w < . The condition (2.5) means that the dif- 
ference between w and 1 is of greater order than R™'”* as R — —; that is, it excludes a 


small range of speeds near the sonic speed. In the range 0 < w < w,, , with the above 
small neighbourhood of w = 1 again excluded, there are two other branches given by 


V = gow) = iE (1 — w’)(1 — au’) "41 + o(1)}, (2.7) 
and 
P R : : 1/2 
V = g,(w) = 2 (l-—w)(1 —- au’) {1 + o(1)}. 


In the range 1 < w < 6°”, there is only one branch V =g ,(w) which is actually the 


continuation of V = g;(w). 
(vii) The curves C, and C, intersect only at w = 0 and w = 6°”. For w < 1, 


f(w) ~ 2w + 2(1 — B)w*® + 2(1 — B)w® + 21 — Bw’ + --- (2.8) 
and 
g:(w) ~ 2w + {2(1 — 8) + 8R"}w" 
+ {2(1 — g) + 8(5 — 486)R™ + 128R°*}w* + {2(1 — 8B) (2.9) 
+ 8(14 — 168 + 36°)R™ + 64(17 — 126)R™ + 2688R™*}w’ + ---. 


Let C, be the curve on which d(V — f(w))/dw = 0 in (1.41); that is, the curve 
given by 
7 a .2 vm 52 2 
a mal - w )(1 Bw’) ‘ oe (2.10) 
(1 — Bw*)* — 4w*R"{1 + (1 — 38)w* + Bu*} 





Then in the range 0 < w < 8”, under the condition (2.5), C; lies between C, and 
V = g,(w). This fact is of use in connection with the approach to the limiting form 3.1. 

2.2. The solution curves may be shown to have the following properties:— 

(i) The point (8°, 0) is a saddle point, the solution curves passing through it 
being tangent to the two branches of C, passing through this point, see 2.1, (iii). 

(ii) For large w there is a family of curves asymptotic to C, , and another family 
asymptotic to V = —43Rw. 

(iii) There is an infinite number of solution curves approaching w = 0 for large 
negative V. These may be shown to have the form 


V = —4}Rw'+A—(4R— 2)w+ OW’), (2.11) 


where A is an arbitrary constant. 
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(iv) By a method similar to that used in (1.3), (or by substitution of a power 
series), the asymptotic form of the solution curves through (0, 0), in the neighbourhood 
of this point, is found to be 


V ~ 2w + {2(1 — 8) + 8R“"}w* + {211 — B) + 8(5 — 48)R" 
+ 128R~*}w* + {2(1 — 8) + 8(14 — 168+ 36°)R"' = (2.12) 
+ 64(20 —158)R~* + 3456R™*}w’ + --- + C’e™’ 


so that by comparison with (2.8) and (2.9), all these solution curves lie above f(w) 
and g,(w) for w small enough, as the coefficient of w’ in (2.12) is greater than the co- 
efficient of w’ in (2.9) by an amount 192(1 — 8)R~? + 768R™. 

The above properties, and a knowledge of the regions of positive and negative curva- 
ture and positive and negative slope, enable the solution curves to be sketched, as in 
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Fig. 2. The w-V plane for constant pu. 


C; is curve of zero slope. 
Cz, 91 » 92 5 gs are the curves of inflexion. 
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Fig. 2. It is seen that the solution curves passing through the origin may be divided 
into three main classes:— 

(I) those which pass through V = g,(w), cross the line w = 8 
asymptotic to C, , for large w; 

(II) Those which pass through V = g.(w), bend round to cross the w-axis, and 
are ultimately asymptotic to the negative V-axis (see (2.11)); 

(III) Those which pass through V = g,(w), and then cross the w-axis and become 
asymptotic to the negative V-axis, as for class II. The solutions of the first class, whose 
curves penetrate into the region of negative temperature are discarded as physically 
impossible just as in inviscid theory. However, it will be seen later that all possible 
flows are automatically confined to 0 < w < 6”? when uz is taken proportional to T. 
3. The approach to the limit. 

Now that the general shape of the solution curves in the w-V plane has been estab- 
lished, we proceed to investigate whether any limiting curves are approached as R > ~, 
and if so, how closely these limiting curves are approximated when R is large. A clue is 
given by considering the differential equation (1.41) in the form 


dV = fe . V(1l — w’) — w(1 - aw’}, (3.1) 


~“/? and are ultimately 


dw wV \2 


which shows that, at a fixed point, dV/dw — © as R — ~. So that considerable portions 
of the solution curves may be expected to have very large slope. In fact, it will be shown 
that these portions are very nearly vertical straight lines. The technique used is to find 
bounding curves which lie on either side of the solution curve considered, and which 
approach a limiting curve as R > o. 

3.1. In Fig. 3 consider the solution curve passing through the point P(b, c) on 
V = g,(w) so that 6 and c satisfy equation (2.2) and 


c= rf i~ Oi = av)" {1 + o(1)}, (3.2) 


but assume that 


(1 — b*)"R” = ofl). (3.3) 


After a considerable amount of calculation based on the value of the solution derivative 
on the bounding curves*, the following results may be established. 
For w > b, the solution curve through P lies above the segment PQ’ of the line 


y wid ae+ a (3 i - F< (1 - av’) —~\w—b, O<e<1, (3.4) 


where at Q’ 
w= b+ 4b(1 — bd’), (3.5) 


and below the tangent 


V = A(w) = c+ B ft (1 — b*) -— ae - av) bw — bd), (3.6) 





*For example, on PQ’ the solution derivative is greater than h’,(w) so the solution curve through P 
cannot cross PQ’. 
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Fic. 3. A typical shock solution in the w-V plane. 
for b < w < 1; in particular, it lies below the segment PQ of the tangent, where at Q 
w=b6+4e1 — db(1 — Db’). (3.7) 
For w < 6, the solution curve through P lies above the segment, PM, of the tangent 
through P, and below the segment PT” of the line V = h,(w), where, at 7”, 
2be « — 2f00)\ 
yah, = - ———_ Tite 
= HI RG — Bl — © ec — ef(d) 
This solution curve also lies below the straight-line segment T’M’, where M’ has the 
coordinates 
w= b, = b,{1 — 2b,(1 — Bb)(A — 1)'7(1 — b/)'”’ RR}, V = f(b,), (3.9) 
with A = 2e"*f(b)/f(b,). 
We denote the points on C, where w = b, w = b, , and w = b, by K, L’ and J’ 
respectively, and as (b — b,), (b; — b,) and (b — b,) are all O(R™'*”"(1 — b*)~'”) from 
(3.8) and (3.9), it follows readily from (2.1) that 


M’J’ = O(R'*7(1 — By”). (3.10) 
Also, if w = b; at the point M, then 
b— b, = O(R'7(1 — b*)"'”), (3.11) 


(3.8) 
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so that when the solution curve through P has approached C, to within a distance of 
O(R~’*(1 — b’)~*”) the distance between the solution curve and the vertical line through 
P is only O(R™*”(1 — b)~™”). 

Furthermore, as the solution curve through P cannot cross V = g,(w) it also cannot 
cross C; , (by the remarks following (2.10)), and so, for this solution curve, v — f(w) 
decreases monotonically to zero as w decreases from b, to 0. Thus, for the solution 
curve through P, the vertical difference between this curve and C, is less than M’J’ 
for w between 0 and B, . 

To sum up the progress so far:— if we keep b fixed and let R > ~, the solution curve 
through P approaches C, for 0 < w < b (since M’ — K), and the vertical line KPN, , 
where N, has the coordinates 


tink. ae + Zc ~vy1—-ofe-f®}, O<p<l. (3.12) 


When p = 1, the distance of the curve from the vertical through P lies between N,Q 
and N,Q’, which are both O(1 — 6’), by (3.5) and (3.7). 

3.2. Let the solution curve through w = 8~'” cross the curve V = g,(w) at w = w,, 
then with reference to 2.5, the curves for which b > w,; > 0 belong to class II and those 
for which 0 < b < w; belong to class I, and by using appropriate bounding curves it 
may be shown that 1 — w; = O(1). (Actually, it will be shown in 4.2 that w; — 6” 
as R — .) Accordingly, if we suppose that b > w, for the solution curve already con- 
sidered, it will cross the w-axis at w = d say, where 0 < d < 6°”. Furthermore, it 
may be shown that, if (1 — b?)"'’R~'”* = o(1), then (d — 1) is positive and at least of 
order (1 — b), and hence (d — 1)~'R™'”* = o(1). This means that b and d differ from 
unity to at least the same order in R, and lie on opposite sides of unity. 

It is seen from (3.1) that a solution curve crosses the w-axis with a vertical tangent, 
and from Fig. 2, the solution curve always lies to the left of this tangent. For the solu- 
tion curve through the point P’”’, w = d, V = 0, first consider the part with 1 < w < d, 
V > 0. This lies to the left of the line w = d, and a bounding curve is required which 
lies to the left of the solution curve. It is possible to find an inclined straight line which 
has this property, but this is not sufficient, as we need later to take an integral of V~* 
over the bounding curve, and the result should be finite. We desire a bounding curve 
which has the same singularity at w = d, V = 0 as the solution curve, that is, which 
behaves like (d — w)'” near this point. The first two terms of the Taylor expansion 
of V near w = d, V = 0 give a satisfactory bounding curve, P’’Q”’, given by 

— 1 


2 
V = hw) = {2(1 — Bd*)Rd"}'7(d — w)'? + “ot R(d — w), (3.13) 


where at Q” 
w = a1 ~i@ - Dp}. (3.14) 


For the part of the curve 1 < w < d, V < 0, the solution curve lies to the left of 
the vertical line segment P’K”, and to the right of the line segment P’M”, where 
M”" has the coordinates 


w = d, = d{1 — 2(1 — pd’)(d’ — 1)"'"R-™},, V = fd). (3.15) 
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Let the points K”, J” lie on C, with w = d and w = d, respectively, and let the point 
Nj’ lie on the line w = d, with 


, Rd’ — 1)’ 108d°(1 — Bd’) |'”” 
vu = s nan: US anal 3.16 
a {1 +| R@ 1) | , O<p<il, (3.16) 
so that N{’ has the same value of V as Q’’. Then 
ee De = 
ON Ns! = ( d(d~ — 1), (3.17) 
) 
2d(1 — Bd") 
Mh " = — am ee ee : Cc 
K"M i? — 1)R? » (3.18) 
and it follows that 
M" J" =OR'*@ — 1). (3.19) 
The solution curve, after crossing M’’J’’, then crosses C, at w = d, say, where 
1 < d, < d, , and thereafter lies below it. Thus in the range 1 < w < d, the solution 
curve lies between C, and V = g,(w), and as the difference in height between V = f(w) 
and V = g,(w) is O(R™* | w’® — 1|~*, (from (2.6)), the solution curves lies within 


O(R'?(d? — 1)°°*”) of V = f(w) from w = 1 + O(R-*(d? — 1)*”*) tow = d, . Con- 
tinued backward from this range, the curve crosses the line w = 1 with V < —(1 — 
6)°"R* (from 2.3) and is ultimately asymptotic to w = 0, as given by (2.11). At the 
point NV’, the difference in w between the solution curve and the line w = d is of order 
R’~'(d’ — 1), so that as R > ~, the solution curve approaches C, , between 1 < w < d, 
and the straight line K”P’’N?’ . 

In the range b < w < d, define the points B, B” to have respectively the coordinates 


w = 1, V=ct+ oh (1 — b*)(1 — b)f{e — f(b}, (3.20) 
2b°c 
and 
w= d, Vee = (1 — b*)(1 — b){e — f(b}. (3.21) 
2b°c : 


Then the solution curve obviously lies above the line Q’Q’’, and below the lines QB 
(the continuation of PQ), and BB’’. Thus the maximum V attained is of order (1 — 6)’R. 

The required bounding curves have now been found for the only important class of 
solution curves, and in the next section, this part of the discussion is completed by 
finding a relation between 6 and d, which shows the shock-character of the solution. 
4. The ‘‘diffuse” shock and the main solution. 

It has been shown in the previous section that over some part of the range of w, 
a typical solution curve of class II approaches the inviscid solution curve in the w-V 
plane, as R — o. It will be shown in this section that the steeply ‘humped’ part of a 
typical solution curve of this class approaches a solution appropriate to an ideal shock. 

4.1. We require first an estimate of the actual physical distance between the points 
P and P”. It is sufficient to consider the difference in ¢ between these points, and since 


V = —2dw/dé, 
then 
t-te =2f 2. (4.1) 
P P *i Ff ; 
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From (3.4), (3.5), (3.12), (3.13), and the rest of Sec. 3, 


nwa’ 





dw . dw : is _ 
Ep — Ep <2 | h.(w) + 2 i = kee ‘(w) + 2(d — b) min (Vq.. , Vo), (4.2) 
and 
tp — tp > 2 | va + 2d — 1)[A()]". (4.3) 


The integrals which appear in (4.2) and (4.3) are respectively 


2b*c AL = 2\2 
aT Pn - a fib) log {1 | 4 bURe(1 — &(1 — OY — s0)}, (4.4) 


- Ca ~ log f 1+ ; d"R'?(d? — 1)°?(3(1 — gd’)]| ah (4.5) 
ki d — 1) \ 
and 
2h*e 
RL bye = fib) log {1 1 + 9 ty *R(1 — bd) — BW) —- so}, (4.6) 


so that, for large R, (by what has gone before) they are all 


flog [ROL — By’ u} 
O RA =b S (4.7) 


Also, by their definitions, Vg. , Vg. and h(1) are all O{R(1 — b)*}, so that the terms 
not involving integrals in (4.2) and (4.3) are all 


O{R“"(1 — Bb}. (4.8) 
Hence 
‘i flog [R(1 — b)’ I} 
Ep — fp = wiles | RG — b) es (4.9) 


for large R. 

4.2. With the aid of the above result we are now able to find the asymptotic relation 
between b and d, which in the limit R = ~, will be seen to be equivalent to Prandtl’s 
relation between the speeds on opposite sides of an ideal shock. 

Equation (1.15), (with constant 4), may be put into the form 


“ (a + B)w + ss = +5 4 {ey la uh, (4.10) 
and so 
tere gy f 4 a : pee’ (60 Aw 
in ae \" + aw + 2 at = — plVe- _ V+] (2 — 4) ae 
2c 


_ R = Op ma Ep) 


as 6/w is O(1) between &p and £’/. From (3.2), the first term on the right-hand side is 


asymptotically 


{2R-*(1 — b*)(1 — Bb*)}"”, 
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while from (4.9), the second term is of higher order, so that 


w=d 
| + B)w + $] = {2R7'(1 — b*)\(1 — Bb*)}'7{1 + o(1)}, (4.11) 


w=b 


and as (d — b) is O(1 — b), 


bd = 1+ 0(\R°'*(1 — b)~”). (4.12) 
In terms of actual speeds, this relation becomes 
UU, = a**{1 + O(R''*(a* — u,))7"}, (4.13) 
where 
/ 9 1/2 
a’ = ( - ) oa (4.14) 
7 — 1, 


and as R — o, this becomes the Prandtl relation. 
With reference to the beginning of 3.2, where the solution curves belonging to class 


II were defined as those for which b > w; > 0, we see from (4.12) that as d — 87", 
b— B?{1 + O(R'*”?(1 — b)~’”)}, and hence we have J 
w; = p"*{1 + OR piss — b)~*”*)}, (4.15) 


4.3. Summary of main solution. From sec. 3 we see that in the w-V plane, a typical 
solution curve of class II which passes through P(w, , V,) on the curve V = g,(w), 
approaches the inviscid solution curve in the ranges 0 < w < w, andl < w < w,. 
For large R and fixed w, we may summarise the behaviour as follows:— 

(i) nO < w < w, — « , V > O, where e, = O(R™”), V goes from 0 to f(w) — m, 
m = O(R~”’), and V — f(w) is at most O(R~™”). 

(ii) nl+e<w<w—e,V <0, (w.=w; [1 +O0(R ”)]), where e, = O(R™””’), 


e, = O(R™”’*), V goes from f(1 + €3) + 73 to f(w2) — m2, where m , m = O(R™™”), 


and V — f(w) is at most O(R~”). 
(iii) (a2) nw, -a Swiwte,”, V > O, where ec,” is O(R”), V goes from 

f(w:) — m to O(R*”). 

(b) In w, — «.” < w<w.,V > O, where e,” O(R~”), V goes from O(R*”) 


to 0. 
(c) Inw, + «,” <w< w, — «”,V > 0, Vis O(R ”), being O(R) when p = 0, 
that is, when both w — w, and w, — w are O(1). 
(d) nw, > w> w,- «,V < 0, V goes from 0 to f(w2) — m. 
(iv) Inl +e«,>w>0,V <0, the curve crosses w = 1, with| V| > (1 — 8)’*°R", 
and becomes asymptotic to w = 0, as given by (2.11). 
The part of the solution described in (iii) above approaches an ‘ideal’ shock solution 
for the equation for the conservation of mass, equation (1.39) and equation 


as R — o~, 
(4.12) give, in the limit, 
Pilly = Pollo , (4.16) 
' y¥-1 2 - 
WP = gi — %—— Vv, (4.17) 
p 2 
and 
2 ; 
Ue = : | Qo ; (4.18) 
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(where the fact that §p — —p.. ~0.as R — = is used in (4.16)), and these imply the 
Hugoniot relations; the curve approached is w = w, , f(w:) < V < © andw =w,, 
f(w.) < V < @. 

The actual shock, (as opposed to the limiting ‘ideal’ shock) has a maximum velocity 
gradient of order R, which occurs where the solution curve crosses C, (within O(R™'’) 
of w = 1) and its ‘thickness’ is of order R~'(1 — w,)~* log [R(1 — w,)*], from (4.9). 
The ‘thickness’ has been defined as the distance in the physical plane between P and 
P” in Fig. 3. Of course, there is some degree of choice as to the definition of the shock 
thickness, and this leads to an apparent discrepancy between the above results and 
the usual result for a plane shock [3], which is that the shock thickness is of order 
up ‘(1 — w,)"'((1 — w,) is proportional to the shock strength.) The difficulty is re- 
solved if an arbitrary length is introduced in the plane-shock treatment so as to allow 
the definition of a Reynolds number. Then, for a given shock strength, it is seen that 
the thickness is proportional to R~™' (implied above) only if the velocity gradients at 
the points measured from are O(R), while if the edges of the shock are taken to be at 
points where the velocity gradient has fallen to less than O(R), as in the treatment 
here, the results agree. In practical cases, however, the difference is slight. 

From (1.21) and (4.16) to (4.18), the solution approached in the physical plane is 


given by 


r= tp wl — aw", OS w sm, (4.19) 
=p Bw, ws, (4.20) 
where 
, _— Aap?) (1-897 (28) 
Po _ wt 1 = | , (4.21) 
Po l— wv," 


and the suffix , refers to conditions at the stagnation point at r = © (pp) and pj are 
different because of the change in entropy through the shock.) From (4.21) we see 
that a change of w, corresponds to a change of the boundary conditions at a given point 
of the supersonic region, and it is this that determines which of the infinite number of 
solution curves of class II is selected for given boundary values. This corresponds to a 
diffuser with constant outlet conditions. 

The remainder of the solution curve in the w-V plane, for 0 < w < 1, with V negative, 
approaches in the physical plane the flow starting from a stagnation-point at 


r=r,= ny i p)7 (8) (4.22) 
oPo 

with infinite pressure and density, and proceeding with infinite velocity gradient up 
to w = 1, where it joins on to the flow given by (4.20) (see Fig. 4). The equations would 
break down before this second stagnation-point is reached; nevertheless, this theory 
does not exclude a transition from subsonic flow to supersonic flow without a contrac- 
tion as, in fact, there is heat addition at the second stagnation point. However, the 
velocity gradient in this range is so great, that it is difficult to see how this kind of 
boundary condition could be realised in practice. 
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Fic. 4. Typical w-r curve for a viscous gas. 


4.4. Some numerical results. This theory could be applied to the flow in the diver- 
gent part of a supersonic-subsonic nozzle, neglecting the boundary layer effects. It 
may be useful for low-density tunnels, as it should be more precise than plane-shock 
theory. To give some idea of the magnitudes involved, some rough results calculated 
from 4.1 and 4.3 are given below, in which 6, and 6, are respectively upper and lower 
bounds for (r,; — 172). The estimates of (r; — rz) could be improved, but not without a 
considerable amount of labour. The thicknesses involved are of the expected order of 





magnitude. 

ry To pi — p2 Po 6; 6 
Ta We cm 0A ra gm.cm.~ R micron micron 
0.7 1.43 17.98 288 1.04 1.226 X 107% 6.02 X 10° 6.12 0.93 
0.7 1.43 17.98 288 1.04 1.226 X 10~4 6.02 * 105 18.7 8.0 
0.7 1.43 17.98 288 1.04 1.226 X 10-5 6.02 X 104 402 67 
0.9 1.11 15.86 288 0.234 1.226 x 1073 6.02 X 10° 14.4 2.57 
0.9 1.11 15.86 288 0.234 1.226 X 10-4 6.02 X 10° 116 20.5 
0.9 ] 


11 15.86 288 0.234 1.226 * 1075 6.02 X* 104 899 157 


4.5. It may be shown by arguments similar to those employed above that if P(w, , V,) 
lies on V = g;(w), and (1 — w;)’*R™’” = o(1), that is, for the curves of class III, the 
solution curves cross the w-axis before w = 1 and then approach w = 0 as described 
before. As R — ~~, the solution curve approaches C,; in0 < w < w, (V < 0), and the 
remainder of the curve gives rise to a sort of ‘negative shock’, similar to that for the 
curves of class II in 0 < w < 1, (VV < 0). 

When (1 — w?) is of order R~’’’, (the curves previously excluded) the situation is 
confused, since there is a transition from curves of the ‘shock’ type to those of the 
‘negative shock’ type just described as P moves along V = g,(w) to (w,,, V,,) and then 
back along V = g,(w). It is impossible to tell without detailed calculation, possibly 
involving numerical integration, just where this transition takes place. 

It should be noted that none of the solutions obtained has an infinite discontinuity 
of velocity gradient, even in the limit R = @, and the flow patterns are quite different 
from those obtained for the inviscid fluid. 
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5. The entropy variation. 

In this section we will consider the variation of specific entropy over a typical solution 
curve of the shock type, as an interesting fact arises which, it is believed, was not hitherto 
known. It is found that there is an entropy maximum within the shock, and in the 
limit, as the shock becomes infinitely thin, this maximum does not disappear. Of course, 
the same phenomenon occurs with the plane shock, as is demonstrated. 

We define S to be the specific entropy, and start from the well-known equation 


7 DS 

Dt 

where ® is the viscous dissipation function. 
For this problem, the equation becomes 





= 7 (8 +0V"7} (5.1) 





Ton 2S = 4 du _ udu , (du) dT , 1dr 5 
pu dr 3 A r dr + (4) } + af dr® + r dr \, (5.2) 
and in terms of the non-dimensional quantities from (1.12), (1.13), (1.14), and (1.22), 
| er ore “| arf 20 , dV at 
Rae NM HQeU ta + BBY aw? + aw dws as 


When o = oa, , (5.3) becomes 


ee an = v3 ~ (1 +o "| _ | — a(t +5 ) |} 5.4 
la” et ra+ea/” |~* rata)” lf 64 


wV dV 
R dw 





igi (5.5) 
The use of (5.5) instead of (5.4) is justified in the same way as the basic approximation 
used in 1.4, as differences in w of order R~* only are neglected. 

In Fig. 3 we see, then, that dS/dé is very small over the regions where the solution 
curve lies close to the inviscid curve, becomes of order 1 in the interior of the shock 
and vanishes very close to the point of maximum velocity gradient. In fact, if w,; and 
ws denote respectively the values of w at which the velocity gradient and entropy 
have their maxima, then it may be shown that 


wy = 1 — O(R(1 — w,)”), (5.6) 
Ws = wy — O(R"(1 — w,)~’) (5.7) 
= 1 — O(R"“(1 — w,)~”). (5.8) 


By examining the signs of V and dV/dw, it is seen that dS/di < 0Om0 < w < wz, 
V > 0, and dS/dt > O over the remainder of the curve. This seemingly paradoxical 
result of decreasing entropy is explained by the fact that as the heat conductivity is 
not zero, fluid elements are no longer isolated systems. In the range 0 < w< ws,V >0 
that is, on the subsonic ‘side’ of the shock, heat is continually being conducted backward 
in the sense of decreasing r, but a given fluid element in this range gives out more heat 
at the rear than it takes in from the front, and so has a net loss of heat energy. 


For a perfect gas, with the neglect of a constant, 


S = C, log (p/p’), (5.9) 
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which becomes, in our notation, 














S — S, = C, log {0(wr)?/"-} (5.10) 
where 
Y C. { 4 A 28 \8) fm 
So = i ar log {RT,« "(1 + £6)*}. (5.11) 
So, at the second stagnation point, S = —o, and with this additional information, 
the sketch curve of entropy variation with ¢ (Fig. 5) may be drawn. If we denote by 
S 
tr r 
Fig. 5. Entropy variation for a typical shock solution. 





S, , S, and S.,,. respectively, the entropy when w = w, , w, and ws , we have from (5.10) 


- y 28/ (1-8) 
S, — So = C, log {6,(wn)”’”°"}, 
S, = Ss, = C, log { 62(We12) 28/(1—8) } ’ 
and 
os \28/(1- 
. = So = as log { 65(Ws?'s) is j ° 


In the limit, as R > ~, w, > w;,', ws > 1 from (5.9), and r, , rz and rs tend to a common 
value, so that 


f ] = % ‘aailaaias| 


Snax — S, = C, log ae ———e Of, ? (5.12) 
max r) 5 ‘= Bw; ie 

y ral y 1 az Bw; 2(1+8)/ —p) & \ 

S, — S, = C, log ie ges B ee ee (5.13) 


the last being the usual result for an ideal shock. 
For plane flow, with constant viscosity, we find the solution 


(w — w,)” 3 m (w, — Wy) } si 
i. | D exp (- -— —~ 2, 5.14) 
(w, — w)"* P\ 4 148 


where m = pu and o = 3/4, [4], which is seen to give a shock-type flow, the ‘speeds’ 
Ww, , W, being attained at z = —o and z = = respectively. In this case 
6 dS _ 3m 1— wo (5.15) 


Ra mite gw %- We. 
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This clearly exhibits the entropy maximum which now falls exactly at w = 1, (the 
point of inflexion of the w-x curve) while the relations (5.13) and (5.14) for the entropy 
differences S,,,. — S, , and S, — S, now hold exactly. 

6. Viscosity varying with temperature. 

In this section we consider briefly the case of » varying directly with temperature. 
The results are qualitatively the same for the shock type curves, but the whole flow 
is now automatically confined to the region of positive temperature. 

The equations to be investigated are (1.15) and (1.16) (with C = 1 as before). 
Assume that 


HK = M09, (6.1) 
and hence 

R = R,/8, (6.2) 
where 

R, = er , a constant. (6.3) 


Equation (1.16) now has the form 


( 
med 40 \, 06 {3¢ d aw | 
6—1+A\1+ Ril + ase Rol +B) \o dé (6 — 1) + 8Bw dé 0, (6.4) 


and to find an integrating factor as in 1.4, we must take 
3 Rk (1+ 8) + 40 


7 4° Rol + B) + 48w?’ 
which is not constant. However, for the range of w in question the departure of o from 
the constant value 3/4 is so small that it may be neglected for practical purposes. With 


this value of oc, if we again define 





(6.5) 


2 


— 
R(1 + 8) 


( 
E ¢—1+ H1 + 


(6.6) 


f 40 ht 
ha a PY | al 


then equation (6.4) has the formal solution 


felt * 8) 46 
E A exp J [ {1 nites) ar|. (6.7) 


Physically, @ must be bounded near the stagnation point at @, and if we assume that 
it is integrable in this neighbourhood, then A must be zero for E to remain finite as 


¢— o, and so, as in 1.4, E = 0, and hence 


1 — Bw’ ' 
Ro(1 + 8) + 48w" 
With change of independent variable to w, and neglect of terms of order R>' compared 
with 1, as in 1.4, eq. (1.16) takes the form 


w V dV 
Ry a aaa Ma 





6 = R,(1 + 8B) (6.8) 


=; Vii — w’) — w(l — Bw’) + 28R5'w*V’ (6.9) 
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and, with 





am V6, (6.10) 
wZ dZ is ~ iol me 
R, dw 2 Z(1 — w’) — w(l — bw)’. (6.11) 


The solution curves for (6.11) are sketched in Fig. 6, and it is seen that in the range 


Z 











Fic. 6. Solution curves for zn « T with Z = VO. 


0 < w < 6”, their behaviour is very similar to those discussed before in Fig. 3. In 
fact, the resemblance between (6.11) and (1.41) is striking, and exactly the same tech- 
nique may be used to discuss this equation. It is obvious that for the ‘shock’ type curves, 
similar results will be obtained for Z as were obtained previously for V, and as the 
variation of @ is only of order 1 for the range 0 < w < 8'’”’, we can carry over the same 
qualitative results for V for this new equation. (Only coefficients will be affected, orders 


of magnitude will be unaltered.) It is only near w = 6 ‘”” that the results are different, 


as 6 vanishes here. To investigate the behaviour more closely, we may now turn to the 


(w-V) plane. 
From (6.9), dV/dw is zero when 


2Bw* V* 


Vil — w’) — w(1l — Bu’) + => = ( (6.12) 
2 k 


’ 
0 
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provided that neither V = 0, w = Onor w = 6~””. Equation (6.12) has the two solutions 


; R, 1 . l . Sew" , Jie 
Vv. = 460 2 5 (1 —w)+ EF (1 — w*)? + — (1 — Bw") | \ (6.13) 


and 





= js, | 5 (1 — w) -|} (1 — w)? + Soe (1 — pw u’) | \, (6.14) 


| 


and provided that (1 — w*)’R, > 1, that is | 1 — w |~*R>'’” = o(1), we have 


V. = f(w)(1 + o(1)), 0<w<l, 
os (6.15) 
_ (w — IR, ; - 
= 4eu (1 + o(1)), l<w<B 
and 
eg oe —¥ v')Ro (1 + o(1)), 0<w<il, 
(6.16) 


= f(w)(1 + o(1)), L¢<qess™. 


Thus V. inO < w < 1, and V_in1 < w < 6, lie very close to the curve C, which 
occurred before, with a small neighbourhood of w = 1 of order R~'”’ excluded (this 
small neighbourhood lies inside the neighbourhood of order R~'”* previously excluded 
for the —— unt curves) . There i is no longer an infinity at w = 1, and V, and V_ cross 
w = | at a height of colar R’ 

For tio V of order R, 


wVdV | 28wV*_ f 46w* \ r 
R, dw (1 — Bu’) \! w+ R, V + o(1) (6.17) 


and we see that in 1 < w < 8°”, d’V/dw’ > O above V., and d’V/dw’ < 0 below V, , 


so that we may now proceed to seoteh the solution curves, as in Fig. 7. The solution 


curves which do not bend round to cross the w axis, now go off to V = +o asw—- 
g°'/? — 0, and in fact, all the possible flows starting with 0 < w < 6°” are confined 
to this range; the line w = 6-'” has become a barrier. The speed at which the accelera- 


tion is 8 maximum is now only within O(1) of w = 1, and lies between w = 1 and 


w = 6 **, but the entropy maximum still lies very close to the point of maximum ac- 


celeration. 

When w, < 8” (in the previous notation), an incomplete shock is formed, starting 
from w = 6 '’*(@ = O) with infinite acceleration. This again leads to an impossible 
boundary condition as regards the production of this flow in practice. Apart from this, 
the flow patterns are fundamentally the same as those for constant viscosity. 

The case of the vortex source may be discussed in an approximate manner with 


similar results to those already obtained. 
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APPROXIMATE ANALYSIS OF STRUCTURES 
IN THE PRESENCE OF MODERATELY LARGE CREEP DEFORMATIONS* 


By 
N. J. HOFF 
Polytechnic Institute of Brooklyn 


1. Introduction. When creep strains of the order of magnitude of one or two percent 
develop during the lifetime of a structure, it is often permissible to disregard the primary 
phase of the creep deformations and to base the analysis solely on the secondary or 
steady phase of creep. For the metals used in structures the experimentally established 
secondary creep law is generally given in the form 

én = Koy ’ (1) 
where ¢;, is the tensile strain rate caused by uniaxial tension in direction 1, o,, is the 
corresponding tensile stress, and K and n are constants. When the creep strains are 
large (of the order of magnitude of 0.01), the elastic deformations can often be neglected 
in the calculations, as will be demonstrated by means of an example. Thus the limiting 
state of stress and strain approached as the creep strain becomes large as compared to 
the elastic strain can be determined on the basis of a simple non-linear stress-strain 
rate law. 

It is believed that structural analyses based on the assumptions stated are satisfactory 
for supersonic guided missiles whose surface is heated to high temperatures by the air 
flow. As guided missiles are generally used only for a single flight and not over long periods 
of time like piloted airplanes, their structure can be permitted to undergo large permanent 
deformations. 

2. The elastic analogue. It will now be shown that the stress distribution in a body 
whose deformations are governed by a generalized version of the non-linear creep law 
of Eq. 1 is the same as that in a non-linear perfectly elastic body provided the elastic 
stress-strain law and the boundary conditions are suitably chosen. Following Prager’s 
suggestions for the representation of the stress-strain laws of strain-hardening materials 
[1]t, the uniaxial stress-strain rate law of Eq. 1 is generalized to read 


E= f(J ’ J3)(p(J 2 ’ J;)T + q(J > ’ J3)S’], (2) 
where 
T = 8S” — (2/8)J.I, (2a) 
E is the strain rate tensor, S’ the stress deviation tensor, and J the identity tensor. 
The typical component of the stress deviation tensor is defined as 
84; te (1/3)8.4.5;; ’ (2b) 





*Received June 15, 1953. Part of the analysis described in this article was performed under Contract 
No. AF33(616)-116 sponsored at the Polytechnic Institute of Brooklyn by the Wright Air Development 
Center, Air Research and Development Command, of the U.S. Air Force. The author is grateful to the 
U.S. Air Force for permission to publish this article. He also acknowledges his indebtedness to Professor 
W. Prager of Brown University for his detailed suggestions regarding the formulation of stress-strain laws 
in triaxial states of stress and the statement of the analogues; and to Professor Frederick V. Pohle of the 
Polytechnic Institute of Brooklyn who was kind enough to check the calculations. 

tNumbers in brackets refer to the References. 
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and s;; is the typical component of the stress tensor. The first invariant of the stress 
deviation naturally vanishes: 


J, = 8, = 0. (3a) 

The second and third invariants are defined as 
Jo = (1/2)s8;;83; , (3b) 
J, = (1/3)8) 858i; (3c) 


In Eq. 2 the symbols p(J, , J;) and g(J2 , J3) denote polynomial functions, and f(J2 , Js) 
an arbitrary function of J, and J; . These functions must be determined from empirical 
data to be obtained from creep tests. Eq. 2 is meant for use only when the strain is 
small as compared to unity (say 0.01); under such conditions it implies that the creep 
deformations are inextensional (e;; = 0). 

If a body which follows this stress-strain rate law is subjected to given body forces 
g;(x, t), (where x is understood to represent the three Cartesian coordinates of space) 
and to given surface tractions 7;(x, t) on a portion S, of its surface while the points 
on the remainder S, of the surface are slowly displaced with given velocities V,(x, ?) 
(so slowly that the resulting inertia forces are small as compared to the forces corres- 
ponding to the stresses S and surface tractions 7’), at a generic instant ¢ the stress field 
8;;(z, t) and the velocity field v;(z, ¢) throughout the body must satisfy the following 
equations: 

(4) 


(0s;;/0x;) + 9; = 0, 
dv; /0x;) + (dv;/dx,;) = 2f(J2 , Js) [p(J2 , Js)ti; + Q(J2 , Js)si;]. (5) 


The three equilibrium equations (4) and the six stress-strain rate equations (5), together 


with the boundary conditions 
8;,n; = T, on aes (Ga) 


= VV, on ie (6b) 


7 = 


define the stress and velocity fields in the body. 
The analoguous perfectly elastic body is required to follow the stress-strain law 


E* = f(J2, Js)[p(J2, Js)T + QJo , Js)8’] (7) 


in which the only new symbol, £*, represents the strain tensor whose typical component 
is e;; . The body forces and the surface tractions are kept unchanged but the velocities 
V,(z, t) prescribed on surface S, are replaced by displacements U;(z, t) equal to V;(2, ¢) 
in magnitude and direction. Under these conditions the displacement field u,(x, ¢) and 
the stress field s;;(x, ¢) in the elastic body are defined by Eqs. 4 to 6 if v,(z, t) is replaced 
by u,(x, t) and V,(z, t) by U,(a, t). Consequently to any solution of the elastic problem 
there corresponds a solution of the creep problem and the stress distribution is the same 
in the two solutions. 

Similar results were published earlier for linear visco-elastic materials by Alfrey 
[2] and Tsien [3]. 

3. Pin-jointed truss. As an example for the use of the analogue the stresses will now 
be calculated in the bars of the pin-jointed framework shown in Fig. 1. One end of each 
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~ a mad 


Fia. 1. Pin-Jointed Framework 


bar is attached to a rigid wall. Bar 0-3, which is in compression when the load W is 
applied to joint O, is assumed to be perfectly rigid and braced in the direction perpen- 
dicular to the plane of the truss so as to prevent its buckling. The material of the other 
two bars is subject to creep in accordance with the law 


e = (a/d)", (8) 


where e is the rate of creep in uniaxial tension, ¢ the tensile stress, and \ and n are con- 
stants. Because of the analogue the creep problem can be replaced by a problem in 
non-linear elasticity with the stress-strain law. 


e = (¢/d)", (9) 

where « is the tensile strain. 

It might be mentioned here that the behavior of pin-jointed structures was investi- 
gated by Meacham [4] on the assumption of a linear creep law. 

3a. Load W prescribed. When the load W is prescribed, the stress distribution can 
be calculated from the requirement that the complementary energy stored in the bars 
must be a minimum. If the force in bar 0-2 is designated as X, it follows from the con- 
ditions of equilibrium that the forces in the three bars are 
Fy. = (VW5/2)W — (V10/4)X, Fo-rw = X, Fo-3 = —(1/2)W — (VW2/4)X — (10) 


As the complementary energy per unit volume (U’/V) is defined as 
(U'/V) = [ eds, (11) 
substitution from Eq. 9 and integration yield 


(U'/V) = [A/a + D](o/)". (12) 


With A designating the common cross-sectional area of the two elastic bars, the strain 
energy stored in the system becomes 


U’ = [(n + Dor Af V5[(5/2)W — (0410/4) XI"? + V2 X**7}. (13) 
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In accordance with the complementary energy principle the differential coefficient of 
U’ with respect to X must vanish. After some manipulations this requirement can be 
written in the form 


(X/W) = (VW5/2)[(4/5)" + (V10/4)]"?. (14) 


After (X/W) is computed from Eq. 14 for any given n, the values of Fy_, and Fo-_. can 
be calculated from Eqs. 10. The results of such computations are plotted in Fig. 2. 
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Fig. 2. Forces in Bars of Pin-Jointed Framework 
When 7 = 1, the linear elastic solution is obtained; and when n — o, the solution 


corresponds to the principles of limit analysis. 

3b. Velocity V, prescribed. If the velocity V, of joint 0 is prescribed in the creep 
problem rather than the load W, in the analogous elastic problem the elastic displacement 
U, of joint 0 must be given. The potential V’ of the unknown reaction force RF at 0 is 


—V’ = RU, (15) 
with both R and U, considered positive when directed downward. The forces in the 
bars caused by the unknown reaction R can be calculated from Eqs. 10 if W is replaced 


by R. The complementary energy can be obtained in a similar manner from Eq. 13. 
The compatibility conditions are in this case 


o(U’ + V’)/axX = 0, (16a) 
o(U’ + V’)/oR = 0. (16b) 


As a consequence of the first of these equations Eq. 14 again holds provided W is replaced 
by R. It follows from the second equation that 


(R/A) = (2/V5)[(4/5)'” + (10/4) ]A(U,/2a)'”". (17) 








1954] STRUCTURES IN THE PRESENCE OF CREEP DEFORMATIONS 53 


3c. More accurate solution of the creep problem. The problem just discussed is so 
simple that it could have been solved on the basis of geometric considerations of the 
deformations without recourse to energy methods. For the same reason it can also be 
analyzed in the case when the deformations are governed by a more complex creep law 
than the one represented by Eq. 8. Such a more complete analysis will now be presented 
in order to show that the initial (linear) elastic stress distribution in an actual structure 
is rapidly replaced by a distribution that, for all practical purposes, is identical with 
the stress distribution derived from the (non-linear) elastic analogue. 

The creep law will be assumed as 


e = (1/E)(da/dt) + (/d)", (18) 
where ¢ is time. If e, is the strain rate in bar 0-1, the rate of elongation dAL,/dt in the 
bar is 


dAL,/dt = V5 ae, . (19a) 


Similarly in bar 0-2 s 
dAL,/dt = V2 ae, . (19b) 
The geometric condition imposed upon these deformations can be derived from Fig. 3. 
1 
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Fic. 3. Displacement of Joint O 


When the length of bar 0-1 is increased from L, to L, + AL, the length of bar 0-3 
remains L,; because this bar was assumed io be perfectly rigid. The new position of 
joint 0 can then be found by drawing an arc of a circle from point 1 with a 
radius L, + AL, , a second arc of a circle from point 3 with a radius L; , and by deter- 
mining the point of intersection of the two arcs. When the deformations are small, the 
ares of circles can be replaced by straight line segments perpendicular to the original 
directions of bars 0-1 and 0-3, respectively. This was done in the figure and the vertical 
displacement A of joint 0 was obtained as 


A = (V5/2)AL, . (20a) 


The same vertical displacement must result if the geometric construction is carried out 
for bar 0-2. Hence 


A= V2AL.. (20b) 
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When these two equations are differentiated with respect to time and their right-hand 
members are equated to each other, the following condition is obtained: 


(W5/2) dAL,/dt = V2 dAL,/dt. (21) 
Substitutions from Eqs. 18 and 19 yield 
(5/4) [(1/E)(do,/dt) + (0,/d)"] = (1/E)(do./dt) + (¢2/d)". (22) 
When the load W is prescribed at joint 0, equilibrium requires that 


(W/A) = (2/V/5)o, + (1/V/2)op . (23) 


99 


Solution of this equation for o, , substitution in Eq. 22, and manipulations result in the 
differential equation 


(dx/dt) = K[(1 — ba)” — (cx)"], (24) 


where 


K = [2"(W/Ad)""(E/d)]/[(5/4) + (8/5)'”]. 


After separation of the variables and integration the following solution is obtained: 


b= (hse) | dx/|(1 — bx)" — (cx)"]. (26) 
A numerical example was worked out in which the material of the bars was 52S-H38 
alluminum alloy. At 400°F the material constants can be taken as n = 5, 


\ = 25,000 hr.’”’ lb. per sq. in. 


on the basis of creep tests carried out by Dorn and Tietz [5]. From Table 3.1211(b) in 
ANC-5 [6] Young’s modulus can be estimated as 9 x 10° lb. per sq. in. With these 
values one obtains K = 487 per hour if the load is assumed to be 22,000 lbs. The load 
is applied at ¢ = 0 and dynamic effects are disregarded. Initially the fully (linearly) 
elastic solution must prevail; correspondingly x = 0.56233 when t = 0. The state of 
fully developed creep is reached only as ¢ approaches infinity; then, in the limit, 
x = 0.61206. Evaluation of the integral in Eq. 26 yielded corresponding values of x 
and ¢ which were plotted in Fig. 4. It can be seen from the figure that after the lapse of 
about 70 sec. x is 0.61 which is only one-third of one percent less than the fully developed 
creep value 0.61206; (the difference between 0.61 and 0.61206 is about four percent of 
the difference between the elastic and the fully developed creep values). As the highest 
stress in any bar is about 15,000 lb. per sq. in. (in bar 0-2 at ¢ = 0), the maximum creep 
rate is about 0.0777 per hour. Hence the total creep strain developed in bar 0-2 in the 
first 70 sec. after the load is applied is less than 0.00151. This should be compared with 
the maximum elastic strain in the bar which is 0.00167. 

Thus the conclusion is reached that in this problem in good approximation the state 
of fully developed creep is reached at a time when the creep strain is about equal to 
the elastic strain. The effect of the elastic stresses can therefore be neglected when the 
state of stress and strain is investigated at moderately large creep deformations (at 
creep strains of the order of magnitude of 0.01). 
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Fre. 4. Variation in Non-dimensional Force in Bar O — 1 with Time 


Finally it is of interest to check whether the time interval of 70 sec. is short as com- 
pared to the time needed for the necking and rupture of the most highly loaded bar. 
In an earlier paper [7] the author derived the formula 


t.. = (1/n)/(e_/X)", (27) 


where /., is the time needed for rupture and o> is the engineering stress at the beginning 
of the creep test. It is reasonable to take o, = 14,100 lb. per sq. in. corresponding to 
the fully developed state of creep in bar 0-2. With this value one obtains ¢., = 3.55 hours. 


Thus there is a long period of time during which the fully developed creep solution is 


valid. 
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ORTHOGONAL EDGE POLYNOMIALS IN THE SOLUTION OF BOUNDARY 
VALUE PROBLEMS* 


BY 
G. HORVAY anv F. N. SPIESS! 
General Electric Company, Knolls Atomic Power Laboratory? 


1. Introduction. The method presented below for solving boundary value problems 
is applicable when the problem can be formulated as minimization of some integral. 
One constructs orthogonal boundary polynomials f, appropriate to certain boundary 
conditions, enlarges them by factors g,—obtained from Euler-Lagrange variational 
equations—to product functions 


Grn = Ind n (la) 


defined over the entire domain, with boundary values 


Grn = x ( Ib) 


and expands the prescribed boundary value—say, ®,—of the unknown function @ into 
the edge polynomials 


>, = bas Cle . (2) 
Then, 
®~ DY crn = Do CnGah n (3) 


constitutes an approximate solution of the problem. An interesting aspect of the method 
is that once the general expression for the polynomials f, has been found, any eigen- 
value and its associated eigenfunction can be determined without prior determination 
of the lower modes. 

Another interesting aspect is that while the functions f, are—by definition—pre- 
cisely orthogonal, with respect to some suitable weight function p, 


Gat) =f ofan ds = bm (4) 


boundary 


the derivatives df,,/ds, Of,,/ds (if s denotes the coordinate along the boundary), while 
not orthogonal, can be regarded—in the problems to be considered—as being approxi- 
mately orthogonal, in analogy to the precise orthogonality of the derivatives 0®,/0s, 
d®,,/ds of the exact eigenfunctions. 

In “first approximation” our approach thus neglects the derivative coupling terms 
altogether and yields the product eigenfunctions (la) discussed above. In “second 
approximation” products of df,/ds, df,,/ds are retained in the variational integral if 


*Received March 25, 1953. 
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m = norn + 1 (orm = n, n + 2 when due to symmetry conditions df,/ds, df,,/ds 
turn out to be precisely orthogonal*), and only higher order coupling terms are neglected. 
In such a case the approximate solution—call it now ¥,—will appear as a linear com- 
bination 


Wn —_ Sa 7 ae + Gant n + a J n+1 (5a) 
of product eigenfunctions, with boundary value 
Yr =Sfnr- 


While the polynomials f, of the second approximation are the same as those of the 
first approximation, the G, functions, again determined from a variational equation, 
will differ from the g, functions. Higher approximations are constructed similarly. 

Functional developments of type (3) where f,(x, y) is some assumed function, 
g(x) is a function determined from an Euler variational equation, have been employed 
previously by Kantorovitch* and Poritsky’. The novel feature of the present approach 
is the systematic development of orthonormal sets f, . The determination of these sets 
greatly simplifies subsequent calculations as one is allowed (in first approximation) 
to regard the various modes involved as uncoupled. 

The method described above was first used to prove St. Venant’s principle for a 
plane rectangular elastic region subject to selfequilibrating edge tractions.” It was 
found that the problem could be formulated as a biharmonic eigenvalue problem. The 
present paper illustrates use of the ‘‘first approximation” technique for a simpler group 
of problems, those relating to the harmonic equation. Section 6 discusses the wave 
equation and illustrates also the method of the ‘‘second approximation.” 

2. The potential is an isosceles triangle. In this section we shall solve the equation 


Vt = 0 (6) 


in the triangular region shown in Fig. 1, subject to the following boundary conditions. 








y 
A 
+h 
= 
ea 2a > xX 
b 
I 
=h 
Fia. 1. 


3As is the case for instance when f, is an even function, f,+: an odd function. 

‘Mathematical theory of elasticity, by I. S. Sokolnikoff, McGraw Hill, 1946, p. 315. 

5Reduction of the Solution of Certain Partial Differential Equations to Ordinary Differential Equations, 
by H. Poritsky, Proc. Fifth Int. Cong. of Appl. Mech. 1938. 

6The end problem of rectangular strips, by G. Horvay, J.A.M. 1953, p. 87. See also discussion of the 


paper in the issue of Sept. 1953. 
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Along the boundaries II, III where 


y= +y = +A(1 — z/D), (7a) 
we prescribe 
®,, = 9, = 0. (7b) 
Along the boundary I, where 
z= 0, —-h<y< +h, (8a) 
we prescribe 
®, = arbitrary function. (8b) 


The integral to be minimized is 


b ptud 
r= | | ((e8/a2)* + (08/ay)"I dy ar. (9) 
The lowest degree polynomial f(y) which satisfies 
~+ub | ath 
f(x, +y,) = 0, (f°(0, y)) = | f(x, y) dy = f°, y) dy = 1 (10) 
J —yb z=0 J —h 


is 


15 )"( vs) 
> = (—— — 4). 11 
f (13 Ys (11) 
The higher degree polynomials, obtained by orthonormalization with respect to the 


lower ones, are 


fi: = (105/16h)'?(y/y.) [1 — y?/yel, 


f. = (3/4)'7[1 — 7y’/vilfo (12) 
fs = (11/4)""[1 — 3y°/yi]f, , 
fs = (91/128)'"[1 — 18(y/y)? + 33(y/yo)*Tfo , - 


Introducing the notation 


of(z, y) ot, af(x, y) = f,, dg =’, [ ff, dy = (ff,) (13) 
Ox Oy da J ys 
and writing 
®m teal te Ain = + An gn(X) fala, y); (14) 


we write Eq. (9) in the form: 


a wees 


I= > A.A, | (fofu}n : {gfulm + fofe tg'f}n > fof. + 9’f}m]dydx. (15) 


“—yb 
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To first approximation this is 


ab 
I= DA; | fg” + Uff.d99' + fz + fg"), dx. (16) 


“0 


To minimize (16), the Euler-Lagrange equations 
1 eo? ae 9 wd 
= (Pg sn — Uf fee + 2fs + fig}, = 0 (17) 


must be satisfied. It is readily seen that 
(fr) = y/h, {(ffi2 + 2+ fi}, =c,/hy , (18) 


where c, is a constant which depends on the order n of the polynomial. The differential] 
equation (17) thus assumes the form 


[1 — x/b)gi]’ — c,g,/h'(1 — 2/b) = 0 (19) 
and has the solution 
g, = (1 — x/b)”, v, = b(c,)'"*/h. (20) 


In particular for n = 0 we obtain 
(f?)=1-% (fofozz + 2foz + fo) = (3 + : VAS - =) (21) 
Jo) = bp? — Jodoze TF “Jor T Jou = lo TOR b, < 


= 3 (7) v3" aad 


and 


In first approximation » is the zeroth eigenvalue, g.f, the zeroth eigenfunction. For 
the special case b/h = 10 there results 


y, = 15.9. (23a) 
This compares with the well-known precise solution 


T 


— - = 15.8 (23k 
2 arccot b/h a8 3b) 


% = 


of the sector of vertex angle 2a = 2 arccot 10. 
3. The potential problem in a semi-infinite strip. The extremely simple problem 
V’& = 0 in the semi-infinite strip of Fig. 2, with the boundary conditions 
@, = prescribed 
6, =0 (24a,b,¢) 
},, = O, = 0 


provides a simple basis for comparing the polynomial method with the exact solution. 
The same polynomials f, are used here as in Section 2, one must only place 


y=h=1 (25) 


In = ” Pies (26a) 
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y 
‘ 
Ir 
+I 
4 ~~ Xx @ 
=| 
I 
Fig. 2 
of the Euler equations 
gn’ — vig. = 0. (26b) 
The resulting eigenvalues are 
y= (f?", wy = 1.58, », = 3.24, ve = 5.05, --- (27) 
as contrasted with the exact values 
nett = 1.57, 3.14, 4.71, «+ (28) 
For the specific example 
®, = ] = y* (29) 
the exact solution becomes 
ge = 1.173e"**” cos 7 — 0.209e7*"*” cos 22H Bes (30) 
while the polynomial method leads to 
gv, = 1.143e7° (1 — y’) — 0.143e7°°*"(1 — 8y? + Ty’) (31) 
In contrast, a one-term Rayleigh approximation 
ge =e "(1 —y') (32) 
gives 
ee =e — 9). (33) 


It is interesting to make a numerical comparison of the solutions Z, A, R. At the 
point (x, y) = (1, 0) they give, respectively 


ge = 0.242, ga = 0.235, gr = 0.167. (34) 

Function (29) can be considered as the zeroth orthogonal polynomial constructed 
from the set 1, y*, y*®, --- . Clearly, the resulting approximation is poorer than when 
the set 1, y’, y*, --- is used. Since both sets are complete’ there arises the question of 


7By Sz4sz’ Theorem (Courant-Hilbert, Methoden der Mathematischen Physik, I. p. 86) the set 1, y~, , 


Y, » °** is complete, in a finite interval, 0 to h, whenever J (Ax)~! diverges. 
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the proper selection of the best set. The answer is: the best set is obtained when the gap 
between successive powers is the smallest, consistent with boundary, symmetry and 
finiteness requirements for the functions and their derivatives. 

4. The two-region problem. Consideration of a simple two-region potential problem 
(Fig. 3) shows an additional advantage in the use of the polynomial method, namely 
that the often troublesome transcendental boundary equations are replaced by a set 
of equations which are linear in the unknowns. Consider 7’ = 0 in each of two regions, 
A and B, with boundary conditions: 

















y 
© 
a 
C 
= X 
-b 
Fie. 3. 
@ = 0 on the edges y=a,y= —b and r=c (35a) 
®@ = arbitrary on edge z=0 (35b) 
®, oe oy 1a k®, = Pz, ’ Diy = Pay, ’ K@ ayy = Payyy ee (35c) 


on the interface, y = 0. The polynomials f,(y) to be used have a different form in each 
of the two regions A and B: 


fao = (y — a)(Aoo + Any), fao = (y + 5)(Boo + Bary), 
fai = (y— a\(Ai+ A y+ Ay’), fai = (y+ b)(Bio - Buy + Byy’), ri (36) 


The four constants in f) are determined from the normalization condition and three 
interface conditions. The six constants in f, are determined from the two conditions of 
orthonormality, and four interface conditions. Successive polynomials are constructed 
similarly. The appropriate orthonormalization condition is 


aa 


fanf am dy + Sanfam dy = Snm - (37) 
For the particular values of the parameters 
k = 3(2)", a=1, b = 1/2” (38a) 
one finds 
vy = (fh?) = 1.75. (38b) 


This compares with the exact value 


vy = 1.67 (39a) 
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obtained by solving 
k tan vb + tan va = 0. (39b) 


5. The potential problem in a sector. The case of the circular sector brings out a 
novel feature in the orthogonality condition. The geometry is as in Fig. 4. Here the 





integral to be minimized for the k-th approximate eigenfunction 


Xi = g(r) f (0) (40) 


I, = | [(dg,/dr)? + (dg,/r d0)*] dé r dr 


=|] | (ifd® + (fi/n*] dr ar. (41) 


For cases where boundary values along the circular are are specified, we are led to 6 
polynomials of type f, considered in Section 3. 
For the sake of variety we list the f, polynomials appropriate to the conditions of 
vanishing normal edge derivative: 
[db/r dé].. = 0. (42) 


They are 


1/(2a)'”?, 


II 


f 
Ji 


I 


(315/316a)'’*(6/a)[1 — 46’/a’], (43) 


I 


(343/384a)'”*[1 — (30/7)(0/a)* + (15/7)(6/a)*], - 


fa 
When radial boundary values need representation, then an expansion in polynomials 


of r becomes necessary. On denoting 


©) = f 96) dr (44a) 


“0 
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and integrating out over the variable r in (41), 


ata 


I, = | [ftXgi/r) + filrg)) ao (44b) 


ag 


we are lead to the orthonormalization condition 
al 
(99:/T) = | (gegi/r) dr = du: . (45) 


#0 


We list the functions g,(r) for the boundary conditions 


0®,/dr = O%,;;/dr = R(r) = prescribed function of r, (46a) 
@=0 at r= Il, (46b) 
@= 6/r=0 at r= 0. (46c,d) 


Condition (46d), which essentially states that the polynomials g,(r) must contain a 
factor r’, results from the requirement that in the corner, r — 0, of the sector the trans- 
verse variation, d6/r 06, of the function must be limited (otherwise @ would not be 
uniquely determined). 

The functions are 


g, = 60'*r*(1 — 7), 


40'*r?(1 — r)\(4 — 7r), 


—) 
il 


gs = 140'r*(1 — r)(5 — 20r + 187’), 


g, = 840'7r?(1 — r)\(4 — 27r + 54r’ — 33r’), --- 


In either case, whether we are concerned with expansions into g,(r) functions along 
boundaries I, III, or expansions into f,(@) functions along boundary II, the remaining 
procedure is the same—insert the polynomials into the integral (41), carry out the inte- 
gration over the proper variable, and solve the appropriate Euler-Lagrange equation 
for the second (unknown) factor of the eigenfunction.° 

In Fig. 5 we plot successive approximations 


R(r) ~ > aag,(r), a, = (Rg,/r) (48a) 

to the function 
Re=r (48b) 
in terms of the functions (47). Since we now have to construct a g, expansion for a 
function R(r) which does not satisfy the boundary conditions prescribed for g,(r), the 


approximation exhibits a pronounced Gibbs phenomenon near the end, r = 1, of the 
interval. 


8Solution of the general problem, when non-homogeneous boundary conditions are specified along all 
edges, is obtained by superposition. 
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6. The wave equation in a rectangle. The second approximation. In this section we 
shall treat the vibration problem of a membrane 


V’°s = K*o" (49) 





R@:' 


R@z SUI-ryl+2r- Le +11) ZL 
Re t4r2(l-ry(i-3r + 30°) 
g) R@= sr il-ry(l +20) 


wn 
Ra Sre(i-r) YA 





9) 














N 
~~ a 
———e 
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Fia. 5. 
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in a simple rectangle region. The example will illustrate the use of orthogonal polynomials 
for a problem where more than two independent variables are involved; it will also dem- 
onstrate the use of the “‘second approximation.” 

Let the problem be as follows: Solve (49) for (x, y, ¢) in the rectangular region, 
Fig. 6, subject to the initial condition 


(x, y, 0) = prescribed (50) 

and to the following boundary conditions (at all #): 
$,, = Py = 0, (51a) 
©, = 36,,,/ax = 0. (51b) 
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Problem (49) is equivalent to the requirement that one minimize the integral: 


a & pt aa 


I=] | | [(@8/az)* + (6/ay)* — K*(a®/a0)" dy de dt (52) 


subject to (50), (51). We write in “first approximation”’ 

P(x, y, ) ~ pe Crige(t, y, ) = z: Crgi(X) fly) te(d, (53) 
where g,(x) and f;(y) are members of appropriate sets of orthogonal polynomials, the 
functions 7;; will be determined by solving Euler’s equations for the time variable after 
the integration over x and y has been carried out, and the C,; are to be determined from 
the initial conditions. 

















y 
! 
a 
a 
I 
4 —— X 
b 
-a 
Ww 
Fia. 6. 


The f;(y) functions should go to zero at y = + a; thus they are the same polynomials 
as those of Section 2, except that y, and h are now replaced by a. The g,(x) must be 
zero for x = O and have zero derivative for x = b. 

This gives rise to the set 


Go = (15/2b)'’*(x b){1 — bl, (54) 
g, = (273/2b)'”*(a/b)[1 — (61/26)(x/b) + (16/13)(x/b)*], --- 
Now, 
(9G) = | gigi dx = bx, 
j 5a,b) 
Gif)= | fifidy = du, 
and, in first approximation, 
(gigt) = (gt) Sx: ; mab 
eG ems (56a,b) 
(fifi) = Fe dur 
so that the integral (52) to be minimized becomes 
[= pe Thi , (57a) 
k,l 


In =f (ot) + UP )irh — Krill dt. (57) 
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The integrals (57b) lead to the Euler-Lagrange equations 


tii + vite = 0, (58a) 
where 
Kvn = gt’) + (f?’). (58b) 
Assuming that initial conditions require that 
7(0) = 1, tii(0) = 0 (59a) 
it follows that 
Ter = COS ry;t (59b) 


The first two values of (gi’)'” and first six values of (f/’)'” are listed below and com- 
pared with the appropriate true eigenvalues (k + 4)m and (1 + 1)x/2 respectively. 

















k b(gf2)'/2 (k + 1/2)x % error 
0 1.58 1.57 0.6 
1 4.85 4.71 3.0 
1 a(fi?)1/2 (lL + 1)x/2 % error 
0 1.58 1.57 0.6 
l 3.24 3.14 3.2 
2 5.05 4.71 7.2 
3 7.04 6.28 12 

4 9.19 7.85 16 

5 11.51 9.42 22 











One arrives at the second approximation y,, by writing, in accordance with (5a, b), 


YiilZ, y; t) 
= Je-rl'e-1,1f 1 + Jul 'r,1-2f 1-2 + Gel erf 1 + Jule. r+2f +2 + Jerr] ear.ati ’ (60a) 


vil, y, 0) = ge(x)f ily). (60b) 


These functions, orthonormal at t = 0, are suitable for expanding the initial value 
(50) of &. To simplify the notation we shall consider the special case of 


Yos = GAT of + Tosfs + TosS 5) + nT iafs (61) 


and to further simplify the calculations we shall neglect the last term, g:713f; of this 
expression. Insertion of (61) into the variational integral (52) and minimization leads 


to the equations of motion 
(90°) + SP )1Tor + (ftf3)Tos + (fifs)Tos + K°Tsi = 0, 
(fifs)Tor + (go) + (fs°)]Tos + (fsft)Tos + K°Tss = 0, (62) 
(fif8)Tor + (f3h8)Tos +al(gos) + (fs°)]Tos + K°Tss = 0. 
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For a natural mode of vibration 


Ty = T, coe rt, (63a) 


Kv? = (g2) + a? = (5/2b)? + a’, (63b) 


the secular equation system becomes 


21 2 ivn 1/2777 5 al 7 
go fe 9 Il T', + 76 7, =f 
r= 99 = 15 no1/2m . 
— 5} 1°°7, + (% —a \r, _ 4 66°T,, = 0 (64) 


? 6'°T, — 2 66'°T; + (788 _ a), = 0 
This yields the frequency equation 
—a® + 192.5a* + 7507.50” + 56306.25 = 0, (65) 
with solutions 
at = 3.141595’, a; = 6.32440’, as = 11.94288”. (66) 
Assigning now to a’ in (64) successively each of the values (66), and solving the corre- 
sponding systems, we obtain 


a, Ss at Oe 1 : 0.131442 : 0.007686, 


ll 


—0.134007 : 1 : 0.333781, (67) 


~ 
— 
a 
a) 
o 

I 


°T, :°T; :°Ts = 0.035561 : —0.329016 : 1, 


where the superscript of 7 refers to the particular frequency a; with which the amplitude 
triplet ‘7, , ‘7; , ‘7's is associated. Consequently 7,; has the form 


Tas ae hs 'T. cos (v,t + ¥i) + As  ¥ cos (vst + y3) + A; we cos (v5t + y;) (68) 


The six constants A, , --- , ys are to be determined from the six initial conditions 
T,,(0) = 0, Tos(0) = 1, —Tos(0) = 
(0) ) 3 1 (0) D, (69) 
T;,(0) = 7T53(0) + 7;;(0) = 0. 
One obtains 
A, = 0.12920, A, = 0.88545, A, = —0.29654, 1=7;=7s=0. (70) 


This then completely determines the function yo; of (61) (if we neglect, as proposed, 
the 7;, term). Separating out from the expression of yo; the coefficient of cos vo3f one 
is led to the second approximation, ¢,% , of the true second eigenfunction 3 (in the 
approximation that the g,f;7,3; contribution is neglible): 
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30 = Ngo CT shi + “Tals + TiJS0 COS V3ot 
= N(x/b)(1 — x/b)(y/a)(1 — y?/a’)[1.941 + 1.654(y/a)? — 18.719(y/a)*] 


- cos {[(5/2b)? + (6.32/a)"]'”"t/K} (71) 


(N is the normalization factor.) 

Since the expansion of the initial value (zx, y, 0) of ® in the functions (60b) is per- 
fectly straightforward (and hence the present procedure does provide one avenue of 
improved solution), one may ignore, until further study has been made, the question 
as to what is preferable: Use simple second approximation functions of type (61) which 

at ¢ = 0—violate slightly the orthogonality condition 


I| (x, y) Py. (x, y) dy dx = xx bi (72) 
or use the conventional approach involving an infinite secular system,’ and extract 
from this system a set of unwieldy but precisely orthogonal approximate eigenfunctions. 
These precisely orthogonal functions will be successive approximations to the true 
eigenfunctions according to whether we retain in the infinite secular system only terms 
along the principal diagonal, also terms along the adjacent sub- and superdiagonals, 
also terms along the next sub- and superdiagonals, ete. 


*See e.g. Eigenwertaufgaben, by L. Collatz, Akad. Verlagsges., 1949, p. 398. 
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A SUGGESTED MODIFICATION OF NOISE THEORY* 
By JULIAN KEILSON (Massachusetts Institute of Technology) 


Abstract. A class of stationary, equilibrium, Markoff processes is demonstrated all 
of which have the same equilibrium distribution, W,(x), and correlation function 
R(t) = E, exp (—t/ro) differing from each other in the number of zero crossings of the 
system per second. The processes are described by an integral equation characterized 
by a parameter y. As y approaches 1, the integral equation passes over into the Fokker- 
Planck equation 


ow _ ow re] 
*° = Bo Fond + Ox 


(xW). 
Since the number of zero crossings per second of the system becomes infinite as y goes 
to one, the degenerate nature of the Fokker-Planck process is made evident. 

1. Introduction. Any stationary Markoffian motion of a system in one dimension 
may be described by an equation of the form 


OW (a, 2b) . 4 > 

~ Sw itl (x, t) | A(z, x’) dx’ + / W(x’, A(x’, x) dz’, (1) 
where W(x, t)dzx is the probability of finding the system in the interval (2, x + dz) at 
time t. A (x, x')dzx’ is the probability per unit time that the system if at z, will jump to 
the interval (x’, «’ + dx’). Equation (1) describes the manner in which changes in local 


probability density can occur. It is seen that 
d . 
| We, dx =0 (2) 
dt J 


so that probability is conserved. 

The stationary Markoff character of the motion is assured by the formulation in 
terms of a time independent transition function A(x, x’). By stationary Markoff is 
meant that if at ¢ = 0, the system is at x , its subsequent distribution is completely 
described by the second order conditional probability: 


W(x, t) = P,(xo | x; 0). (3) 
For some A(z, x’) the process will be an equilibrium process in the sense that 


lim P.(2» | x; t) = Wo(2x) (4) 


is independent of 2 . 

If A(z, x’) is not localized about x = x’ thestochastic motion described is discontinuous. 
The system jumps about from point to point. As A(z, x’) becomes more localized about 
x = x’, the jumps become on the average smaller, i.e., the motion becomes “more con- 


tinuous’’. 





*Received April 7, 1953. The research in this document was supported jointly by the Army, Navy 
and Air Force under contract with the Massachusetts Institute of Technology. 
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In this paper a class of transition function A,(z, x’) will be studied. It will be shown 
that, as y — 1 and A,(z, x’) becomes increasingly localized about « = 2’, the integral 
equation (1) passes over into the Fokker-Planck equation 


a W rs) , 
—,- (5) 


ow > Oo, 
E, 37 + ax (xW). 


To 4 _ 5 
Ot Ox 
The nature of the motion described by the Fokker-Planck equation as indicated 
by this limiting process will be discussed, and the physical implications of this study 
examined. 
2. Properties of the integral equation. Consider the process described by the transi- 


tion function 


: ] ; 
Alz, x’) = T ~ (8 m)'/? exp [—A(2’ — yzx)"]. (6) 

Here, 7',.; , the mean free time of the system, is independent of 2, i.e. 
| A(z, 2’) de’ = 1/T ay (7) 


is the mean number of transitions of the system per unit time; 6 describes the dispersion 
occurring at each transition. The larger 8, the more “continuous” is the motion; y is a 
parameter describing the relaxation or damping of the system, i.e. if (x) is the mean 
value of z, 


> 


(7) = | xW (a, t) dx (8) 


d 


then 
d(x) —(z) 
(9) 


dt T,,/4—¥y) 


as is readily deduced from Eq. (1). The decay time of system is, then 
7=T,,/(1 — y) (10) 
The closer y is to one, the smaller is the damping and the longer the decay time. 


7 is also the correlation time of the system, for 


RG) = / W (2%) XoP’s(Xo | 2; Hx dxo dz. 


Since from (9) (x) = 2» exp (—t/r), we have 
R(t) = / W .(ao)xo exp (—t/r) dxo , 
i.e. 
R(t) = E exp (-—?/r), (11) 


where 


Su / W(2)2? de. (12) 
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Equation (1) with A(z, x’) given by (6) may be solved in the following way. Use 


is made of the expansion 


exp [—(x — vy)?/(1 — °)] = {[r(1 — 7°)” exp [+(%* — 2*)/2]}. 
De Vax) aly), 
where y,(z) is the set of orthonormal Hermite functions 
1 1/2 * . 
y,(x) = (<tr) exp (x”/2)(—d/dz)" exp (—2’). 
Then 
2 @ 
A(z’, xz) = “a exp | -< (x? — 2”) | > y'v.lax)y,(ax’), 
mf 0 
where 
a= [a(1 = 7’)]'”. 
In terms of 
a 
¢(x, t) = exp E 2 |we, t) 
Eq. (1) becomes 


@_ > </> Natal ae! , 
alee 20 + tT) % Vn(ax) Prlax’)y"o(x’, t) dx’. 
If we expand ¢ in terms of the orthonormal set a’y,(ax) 


¢(z, t) = p> a,a'’* (ax), 


our equation separates into 





da, _ __G , YQ 
dt i <7 
If, att = 0, 
W(x, 0) = hx — %) = Y avalaz) Yolo), 
then 


a,(0) = a’” exp (= 13) ¥.(oa), 


and 


2 


P.(2o | 2; t) = exp E (x5 — “| > fay,(ax)p,(ax)) exp [—t(1 — y")/Tays]}- 


(13) 


(14) 


(15) 


(16) 


(17) 
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We see that 
lim P(2 | zt; t) 


t—@ 


W(x) = 2 a exp (—a’x’) 


m*18(1 — y°)]' exp [—B(1 — y')2"] (18) 


and 
E = (x’) = [28(11 — 7y’)J"'. (19) 
3. Passage to the Fokker-Planck process. The Fokker-Planck equation, (5), con- 
tains two parameters, the correlation time 7, and the equilibrium mean square, FE, = (zx’). 
Our transition function A(z, x’) contains three variables 7, E, y. Let us maintain 


the values r = 7) , and E = E, and permit y to pass through a set of values approaching 
one. To do so we adjust T,,, and 6 to the value of y through the equations 


Ted = T(1 = 7) (20) 


and 





= ] na ie ] : = 
B - 2E,(1 sg 7’) a 4E (1 7 y) (21) 


These parameter values define a set of processes A, all of which will have the same 
equilibrium distribution and the same correlation function as the corresponding Fokker- 


Planck process. 
Indeed, the corresponding Fokker-Planck process is just the limit of the process 


A,asy — 1. 
This may be seen in two ways. First one may pass to the limit y = 1 in the condi- 


tional probability function. 
Since for the process A, 


7 1/2 
E, . l-y n 
a=\|>) ; lim ws 


7 
mf To 





we have 


lim P2,(%» | x; t) = exp E (% — 2) | D {ay,(ax) p,(ax») exp (—nt/r) 
7-1 _ 0 
—[xz — 2% exp iol (29 


2E,[1 ba exp (—2t/ro) | 








= {2E,r[1 — exp (—2t/7,)]}~'”? exp { 


and this is indeed the second order conditional probability of our equation (5). 
It can also be seen directly that the integral equation passes over into the Fokker- 


Planck equation. 
The integral equation can always be rewritten formally* as a partial differential 


equation of infinite order 


We, wa, 2 ‘ 
t= 2 gg (Antz) WC, 0), (23) 





*See p. 246 Keilson and Storer, Q. Appl. Math. 10, 243-253 (1952). 
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where 
Ala « Se / (x — x')"A(a, 2’) de’. (24) 


In the limit y — 1, A,(x) and A,(z) do not vanish, but all higher moments do vanish. 
It is readily found that 





A,(z) = 2/t 
and 
_(—-y’ ,d-yE .E 
A 2(2) a ZT wat + QT ns c To ; 
The higher moments 
: 1 F 4 1 1/2 ‘ = , 
A,(t) = = / {((l — ya + (ya — 2’)} T., (4) exp [—B(x’ — yx)*] dx 


contain terms 
(1 — ¥)*2"(yz — 2’)’), 


where either m > 2 or p > 2. A simple examination reveals that all such terms go to 
zero as y approaches one. 

4. Zero crossings of the system A,. It is easily seen that the number of times per 
second that the system with its motion characterized by A(z, x’) will cross zero is given by 


0 


j(0) = | Wola) dx [ A(z, 2’) de’ 


_ a2\l/2 pO po 
-G—y' f [exp [—2" — o? — 2yay] de dy 


TT ms J-o /0 


-* =a I, exp [—S*(1 + y) — T°(1 — y)] ds dt, 


T 


A 


eee aa ae 
- 














e) S 


Fic. 1 
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where A is the shaded area shown in Fig. 1, i.e. 

l - 1-—y)\'” = 

0) = —— tan™ {(=2) \ (25 

oe i+7 


For the process A, , 


~3 (3 —_ x" oR 
= - — ts — : 2K 
7+(0) Sail me tan i ‘ (26) 


As y approaches one, 7,(0) becomes infinite. This will be true not only for zero but 
for all x. The implication is plain. The Fokker-Planck process is a degenerate process 
in which the one sided current density of the system is infinite. A Fokker-Planck model 
for the velocity motion of a colloid particle would describe an infinite number of changes 
of direction of the particle per unit time. Such a model used to describe voltage fluctua- 
tions would imply an infinite number of polarity reversals per second. Since a process 
Ay will afford the same correlation function and equilibrium distribution, and finite 
polarity reversal frequency, it is suggested that such a model may better describe noise, 
and that the number of zero crossings be regarded as an independent macroscopic 
physical quantity on an equal footing with 7, , Ep . 

The author thanks Dr, Franz Stumpers of Phillips Eindhoven and Dr. Edwin Akuto- 
wicz for. their interest and encouragement. 


EVALUATION OF CONSTANTS IN CONFORMAL REPRESENTATION* 
By SAMUEL I. PLOTNICK anp THOMAS C. BENTON (Pennsylvania State University 


In using the Schwarz-Christoffel transformation [1], 
de= K [] (¢ - ¢,)*”” "ag = Kf) at 
t=1 


whereby the upper half ¢-plane is mapped into a simple connected polygon, the evaluation 
of the unknown constant K (if complex K = ce"”, c, d real), is oftentimes tedious. We 
shall show a simple method of evaluating the unknown constant K by examples, proving 
first a 

THEOREM: By the Schwarz-Christoffel transformation if ¢; in the ¢-plane corresponds 


to two points P; , Q; in the z-plane and ¢ = ¢; is a simple pole of f(¢), then 


, _ dist (P; , Q:) 


~ miR(e = £;) 


R, denoting residue and dist (P; , Q;), denoting the distance between the two points P; and Q; . 


*Received May 8, 1953. 
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Proof: Since 


dz = Kf(g) dt 


Qi 
[ dz = dist (P; , Q;) 


’ P; 


=Lim | Kf(ts + b¢")ise* do 


6-0 0 


= K Lim | f(t + de'*)ide™ de, 


6-0 #0 
% 7 ¢ 16 
where { — {, = 6¢ 
Since f has a simple pole at ¢; the Laurent expansion [2], is 
R 
a - g 
f(O) rf, g() 


with g(f) analytic near ¢ = ¢; . Now, 


fbi + be) = <5 + ote + Be). 


Therefore, 


dé + i6 [ Law” ao 


0 k=0 


Lim | f(t + be*)ide"’ do = Lim | i [ 


8-0 +0 6-0 0 


imR 


whence 
_ dist (P; , Q:) 
— wtR(F = f,)° 


Consider the transformation described by Milne-Thomson [3], whereby we map the 
infinite strip on the upper half ¢-plane. The Schwarz-Christoffel transformation gives 


ant i (fs 
ar = XS or z=K , +b. 
L = 0 for z = 0 corresponds to ¢ = 1. 
By the theorem, 
- _ dist (P;,Q) _ _a@_ _a 
~ miR( = 1) wil) 


In this case comparison of real and imaginary parts with infinities is avoided. 
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A UNIQUE LAW FOR IDEAL INCOMPRESSIBLE FLOW WITH 
PRESERVED PATTERN OF FINITE SEPARATION* 


By H. 8. TAN (Cornell University) 


Prandtl, in treating the rolled up separation surface for a flow around corner (Fig. 
1 (a); Ref. 1), and von Karman, in treating the wake formation for a flow normal to a 














(a) (b) 





Fic. 1. (a) Prandtl’s case of flow around a corner. 
(b) von Karman’s case of flow past a flat plate. 


flat plate (Fig. 1(b); Ref. 2), both arrived at a law of fluid motion with stationary flow 
pattern which can be reduced to the form 
: Uo 

wy 1 (dot/U>)’ 
where we denote by U(#) the velocity at r = 1 for the first case, and at infinity for the 
second case, by a the acceleration at the same point; ¢ is the time, and the subscript 0 
indicates the initial condition. 

It can readily be shown that this law of motion is in fact at once general and unique. 
The law is general in the sense that it applies to any body shape with well defined separa- 
tion points, provided there exists a solution. The law is unique because no other fluid 


motion except this one will produce such a preserved flow pattern. 
To prove these statements, let us observe that a satisfactory non-stationary velocity 


potential must be of the form 
P(x, y, t) = Ulid(a, y) (1) 


in order that the flow pattern may not alter with time. This leads directly to the following 
expression for the pressure: 


“of toflO +9 
p= 9 4 Ly (2 T \ay : (2) 


*Received May 8, 1953. This research was supported by the United States Air Force under Contract 
No. AF33(038)-21406, monitored by the Office of Scientific Research, Air Research and Development 


Command. 
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Now, ¢ is determined by the differential equation 
a a 
co 45 (3) 


dx” * dy’ 


together with the boundary conditions 














plate flow corner flow 
i) at infinity ¢,=1 ¢ =0 
ii) on solid boundary d, = 0 dn os ea () = 
iii) on free boundary p =0 = 0 





The only boundary condition that involves U, and hence ¢, is (itz). Since determina- 
tion of @ cannot depend on ¢, (ii7) must be reducible to a product form, i.e. 


L(U)M(¢) = 0. (4) 


But this is evidently the case when and only when U satisfies the following differential 


equation 
dU os 
es ae 1) q 5 
FT +kU’, k>0 (5) 
by virtue of (2). 
Integration of (5) leads immediately to the desired law of fluid motion: 
Re «ie U,>0 (6) 
1 + kU,t’ ; 
where U, is the initial velocity at r = 1 for Prandtl flow, and at infinity for Karman flow. 
Introducing the initial acceleration at corresponding points a, , we see that (6) further 
reduces to 
” Us 
(6 = 7 
= T= (at/T) © 
which result has been plotted in Fig. 2. 

It may be observed that a, is uniquely determined in terms of U, and a reference 
length 1 (where ¢(1) = 1 in Prandtl flow, and is the plate width in Karman flow) by 
certain kinematic considerations. Indeed, for the case of Karman’s plate flow, the 
condition that the wake must form a closed region leads to the following value for initial 
acceleration at infinity 

a) = 0.298 Ue (8) 
Whereas for the case of Prandtl’s corner flow the condition that the rolled up vortex 
sheet must coincide with the stream surface leads to the following result 


exp (273'”) —1 U5 (9) 


2{1 + exp (—73'””)] 1° 





a = 
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From Fig. 2 it may be noticed that although it is entirely possible to produce a 
preserved pattern of finite separation both by an accelerated and a decelerated flow 
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Fic. 2. Variation of flow speed U with time ¢, according to Eq. (7), for accelerated (above) and 
decelerated (below) flow. 




















[+ 


field, yet to maintain such a flow pattern indefinitely, an accelerated field is never 
adequate; the flow then must be a retarded one. 

This analysis definitely rules out the possibility of preserving a finite wake in a 
stationary flow field, or in any flow field that does not exactly follow the law of motion (7). 
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THE BOUNDARY LAYER ON A QUARTER INFINITE FLAT PLATE* 
By L. TRILLING (Massachusetts Institute of Technology) 


This note discusses the incompressible boundary layer on the surface of a quarter 
infinite flat plate in the absence of a pressure gradient, generalizing the classical two 
dimensional Blasius solution [1] and Sears’ extension to an arbitrarily yawed plate [2]. 
It shows that the flow retains free stream direction and Blasius profile at all points of 
the plate, and that the projection of the constant velocity surfaces on planes parallel 





*Received May 21, 1953. 
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to the plate are curves which become asymptotically parallel to the plate edges far 
from the lead corner so that the Blasius and Sears solutions are asymptotic cases of 
the solution given below. 

Let the plate occupy the first quadrant of the X, Y plane with the lead corner at 
the origin; let the free stream velocity vector U.. be directed along the line X — Y = 0 
(the results obtained in this case can easily be generalized to an arbitrary angle of ap- 
proach by suitable stretching of the Y coordinate). 

The equations satisfied by the flow in the resulting boundary layer are 


Ux + Vy + Wz = 0, (1a) 
UUs + VUy + WUz = Wer, (1b) 
UVxy + VWVy + WVz = rVzz, (1e) 
with the boundary conditions 
U=W=V=0 on the plate, (2a) 
U=V=22"°U. as Z— ~, (2b) 


If one seeks solutions of the form suggested by Sears and satisfying the boundary 
condition (2b), namely with 


U(X, Y, Z) = V(X, Y, Z), (3) 

Eqs. (1) become 
Ux + Uy + Wz = 0, (4a) 
U(Ux + Uy) + WUz = WW zz . (4b) 


It is convenient to introduce the Blasius parabolic coordinates 
a= Z(Xv/U.)'"*, y= 24(Y¥v/U.)’, (5a) 
and the dependent variables 


9/377 1°WZ 


u(x, y) = > w(x, y) = (5b) 
The equations of motion then become, with r? = x’ + y’, 
—rw, + w+ 5 (z*u, + y*u,) = 0, (6a) 
wru, — (*u, + y*u,) = r°u,, , (6b) 
with the boundary conditions 
u(0) = “2 =o, (7a) 
lim u(z, y) = 1. (7b) 


z.v-2 
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Note that with y = 0, 8/dy = 0, r = 2, the flow becomes a Blasius flow since one has 


l ; _ 
mw,—-w+s5ru, = 0, (Sa) 
u 3 2 *] 
WLU, — SLU, = Tz; - (Sb) 
Differentiating (8b) and substituting for w, w, from (8a, b), one obtains 
(u’’ /u’)’ + u/2 = 0 (Sc) 


which is equivalent to the Blasius equation. 
To reduce system (6a, b) to system (8a, b) in one variable, one must find a parameter 


s(x, ) which satisfies the conditions 


0 fa) 0 0 ‘ 
s— =F =7- y—, (9a) 
Os 0} Ox “ Oy 
fa) ; 0 ; 3 ra] 
[—= tS eS. (9b) 
Os Ox “ Oy 


In polar coordinates, the condition (9a) is satisfied if 


s = rX(6), (10a) 
while (9b) is satisfied if s = rd(@) is a solution of 
as. +ys, =s (10b) 
or, using (10a), 
\(cos* 6 + sin* 6) — d’ sin 6 cos 6 = 2d’. 10c) 
In terms of the new variables 
a = cos 48, oe i (lla) 
Eq. (10c) becomes 
u(3 + a) + 21 — a )p’ = 4p. (11b) 
A particular solution of this equation is 
Mo = (1 — a)/4. (12) 
The general solution can be obtained by substituting 
up = F(a)p(a) (13a) 
into (11b). The function F(a) then satisfies the separable equation 
(13b) 


2(1 + a)F’ = F(F — 1), 


so that one has: 
F =([l1—c(1+a)'"J", (13¢) 
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where C is an arbitrary constant. When C = 0, one finds the solution u.(a). Substituting 
(13c) into (13a) and (lla) and then (10a) 


s = (1/2)rsin 26[1 — a cos 26)”. (14) 


Since s(r, @) satisfies equations (10a, b), u(s), w(s) satisfy the Blasius Equation 
(8a, b). Writing s in terms of physical coordinates, we have 


r= Z[(U.(x + y)/vry]'”, sin 6 = (x/x + y)'”, cos 6 = (y/x + y)'”, (15a) 
s = 2Z7(U./v)'"[(a@ + y) — a(x — y)]'”. (15b) 


The constant velocity surfaces are the surfaces s = const. since u = Bl(s). It remains 
to determine the constant | a | in such a way that as X/Y — o, s becomes proportional 
to Z(u../vy)'’* so that one obtains the Sears yawed plate solution. This is the case if 
a=-l(t@>y)ja=1(@ < y). 

[t follows that the present solution extends Sears’ result up to the immediate corner 
of the plate with no change. The constant shear and boundary layer thickness lines are 
parallel to the axes up to the axis of symmetry. This indicates the existence of a narrow 
region there where the cross flow terms in the velocity Laplacian are not negligible and 
the boundary layer equations do not hold. 
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NOTE ON MAXIMUM SHOCK DEFLECTION* 
By GARRETT BIRKHOFF anp JOHN W. WALSH (Harvard University) 


The angle \ = A(M) of maximum shock deflection, for a given Mach number M 
of flow, is of interest in various applications. It gives the critical angle for attached shocks 
past a wedge [2, p. 53], and that for jetless wedge collapse [3]. We give here a new simple 
means of determining \(/), for ideal gases, which seem simpler than the usual one 
[1, $122]. 

We follow the notation of [2]. By formulas (4.3) and (3.3) of this reference, 


hess a 


qa* 2 
y¥+1qn 


€ 
a= and j= 


» 
= 1 
Jin i. YT (1) 


in a polytropic gas, with p + po = Ap’. 

Now suppose we are given the velocity gq, , relative to J, of one impinging stream, 
and the shock angle 8 between the impinging stream and the shock front, as in Fig. 1. 
Of course, 8 is not known a priori; we shall seek, by variation of 8, that value of 8 which 
maximizes the deflection angle 6, and thus obtain the desired maximum deflection 
Snax = X. 


We suppose also that we know p, p in the impinging stream. Then the normal shock, 
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84 NOTES [Vol. XII, No. 1 











Fig. 1. Stream deflection by a shock. 


associated as usual with the oblique shock by considering those velocity components 


that are normal to the shock front, has on the supersonic (approaching) side a velocity of 


din = Q, Sin B. (2) 
Substituting from (2) and (1), we obtain the important formula 
Oo» a®* FP a. j y-—1 2. - 
: = B “~ = sin 6 2 71 + stig esc” 6 (3) 
Equivalently, since Mi = q;/a’, 
"4 = Asin8 + K esc 8, (4) 
where A and K are given by 
A=(y—-1)/7t+), K = 2/M, (vy + 0), (5) 


for an ideal gas. 
Using trigonometry, we now find that the component of q, parallel to q, is (ef. Fig. 1) 


q cos’ 8 + Jon Sin B. (6) 
The component of g, perpendicular to q, is similarly 
g, cos Bsin B — Qo, Cos B. (6’) 


Taking the ratio, after dividing through the numerator and denominator by q,, and 


using (4), we get 
cos B sin B — cot B(A sin’ B + K) (7) 
: dt oe. F 2.~: Id ; 7 


e 


tan 6 = . —y = 
cos B+ Asin B+ K 


We next rationalize the right side of (7), by the substitution r = cot 8. This gives 
7’ + 1. Substituting in (6’) after dividing through the numerator and denomina- 


a 
esc’ B = 


tor by sin’ 8, we get 
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aaatnel 1 4— Kiet 8). E -B| 
tan 5 = 7] A RED “Ir+or *) (8) 





where 
B= K/(K +1), E=(1—-— A)/(K+)), F=A+K, 
G=K+1. (8’) 


This is equivalent to [2, formula (4.26)]. This is, incidentally, a very convenient formula 
to use in computing shock polars. 


r7O 
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Fic. 2. Maximum stream deflection \, as a function of y and the Mach number M of the entry flow. 


It is especially convenient in computing the limiting angle dns. = . Since tan 6 is an 
increasing function, this is equivalent to maximizing the right-hand side of (8). But 
this maximum may easily be found by setting the derivative equal to zero. The derivative 
is 

E = 2EG?’ as 
F+Gr (F+Gry 





= f* 
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Multiplying through by the denominator, we see that this derivative vanishes if and only 
if r° satisfies 


BG’r* + (2BFG + EG)7’ + (BF’ — EF) = 0. (9) 


Our interest is confined to M > 1, y > 1. Then it may easily be shown, using (8’), 
that the two roots of the equation are real; further, one of these roots is negative and 
has no physical significance since it leads to an imaginary value of r. 

A graph of \(M, y), computed by using formula (9), is shown as Fig. 2. 
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ON THE THEORY OF THE BULGE TEST" 


By E. W. ROSS, JR. anp W. PRAGER (Brown University) 


Summary. It is shown that the use of Tresca’s yield condition and the associated 
flow rule leads to a simple theory for the bulge test for perfectly plastic or strain-hardening 
materials. The basic equations can be integrated in closed form even for finite deflections. 

1. Introduction. The ductility of sheet metal under balanced biaxial tension is 
determined by the bulge test: a circular sheet of uniform thickness is clamped round 
the periphery and subjected to unilateral fluid pressure which causes the sheet to bulge 
plastically. The strain at the pole of the bulge is measured by means of a grid inscribed 
on the originally flat sheet, and the stress at the pole is computed from the applied 
pressure and the measured curvature and thickness of the deformed sheet. The dimen- 
sions of the sheet are chosen so that its flexural stiffness is negligible’; on the other hand, 
the strains cannot be treated as infinitesimal. 

The first consistent theory of the bulge test was given by Hill.* This theory is based 
on the yield condition and flow rule of v. Mises.* It is assumed that the relevant quantities 
can be represented as power series in the ratio between the maximum deflection and the 
radius of the die aperture through which the sheet is made to bulge. Powers higher than 

‘Received Noy. 13, 1953. The results presented in the paper were obtained in the course of research 
sponsored by Watertown Arsenal Laboratory under Contract DA-19-020-ORD-2598. 

2Even for a very thin sheet neglecting the flexural stiffness may not be justified in the neighborhood 
of the edge. Such edge effects are known to be highly localized, however, and may therefore be neglected 
in the discussion of the states of stress and strain in the neighborhood of the pole of the bulge. For a 
discussion of plastic edge effects the reader is referred to a paper by F. K. G. Odqvist (Reissner Anniver- 
sary Volume, J. W. Edwards, Ann Arbor, Mich., 1949, p. 449). 

3R. Hill, Phil. Mag. (7) 41, 1133-1142 (1950). 

4R. v. Mises, Goettinger Nachrichten 1913, 582-592 (1913). 





1954 E. W. ROSS, JR. AND W. PRAGER 87 


the second are neglected in Hill’s analysis which is therefore restricted to moderate 
deflections. 

The present paper contains an alternative theory of the bulge test. This theory is 
based on Tresca’s yield condition® and the associated flow rule’; its equations can be 
solved in closed form without the use of special assumptions concerning the magnitude 
of the deflection. 

2. Yield condition and flow rule. The middle surface of the bulged sheet is a surface 
of revolution, and the principal stresses at a generic point of this surface are directed 
along the parallel circle, the meridian, and the normal. In this order the principal stresses 
will be denoted by o, , ¢,, , and a, ; the initial thickness of the sheet will be denoted by 
h, and the radius of the die aperture by a. For the usual small values of the ratio ho/a, 
the stress o, is much smaller than the stresses o, and o,, . The state of stress at any point 
of the sheet is therefore essentially one of biaxial tension with the principal stresses 
oO and o ; 

In the following, elastic strains will be neglected and the sheet material will be treated 
as incompressible. The principal plastic strain rates will be denoted by e, , €, , and €, , 
and the yield stress in simple tension by oc. 

Tresea’s yield condition specifies that, for plastic flow to occur, the maximum shearing 
stress must have an intensity, o/2, that depends on the state of hardening of the material. 
If the three principal stresses are unequal, the state of flow is supposed to be one of pure 
shear in the plane determined by the largest and smallest principal stresses. For the 
states of biaxial tension considered here, the following basic possibilities must be con- 
sidered: 


a) Ifo. > o, > 0 during plastic flow, then 

Gc. =e¢ and e,. = (Qj, ¢ = —~e, > 0: 
b) if ¢, > ¢, > O during plastic flow, then 

Can = 6 and e, = 0, én = —€, > 0. 


Finally, it is assumed that any combination of the flow mechanisms specified under a) 


and b) is possible when ¢. = o,, . Thus, a third case is added to the list: 
c) if ¢ = ¢, > O during plastic flow, then 
o e = ¢ and e, > 0, ., > ©. €, = —(e. + €n)- 


The flow rule c) is a natural extension of the flow rules a) and b); in fact, it represents 
the only way in which a continuous transition from a) to b) can be achieved. 

For the purpose of comparison, we state the yield condition and flow rule of v. 
Mises for biaxial tension with the principal stresses ¢, , ¢,, and the principal strain rates 


+ Ge 
9 9 9 
Fo — 6.0m + On = OC, 
9) 
€. My = 6 , . ; 
=5 - sign e. = sign (20, — o,). 
‘i. 26m —~ Gc 


5H. Tresca, Mémoires prés. par div. savants, 18, 733-799 (1868). 
‘W. Prager, On the use of singular yield conditions and associated flow rules, J. Appl. Mech. 20, 
317-320 (1953). 
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This yield condition and flow rule has the great advantage of being valid throughout 
the entire range of biaxial tension, thus avoiding the necessity of discussing separate 
cases such as a), b), c) above. On the other hand, each of these three cases has the advan- 
tage that the values of two of the quantities o, , o,, , €. , €m are known outright. This 
fact is responsible for the considerable mathematical simplifications obtained by the 
use of Tresca’s yield condition and the associated flow rule which was first pointed out 
by Koiter.’ 

3. Kinematical considerations. It will be shown in Section 4 that the formation of a 
spherical bulge of uniform thickness is compatible with Tresca’s yield condition and the 
associated flow rule. In this section, the kinematical relations will be derived which 
apply during the formation of such a bulge. 

Figure 1 shows the circular meridian of the middle surface of the bulged sheet: 


ee ee Oe al 
P ie 5 t 








\-o> 





Yo 
Fig. 1. 


O is the center, R the radius, and A a generic point of the meridian; P is the pole, and 
E represents the edge. The angle POE will be denoted by @ and the angle POA by ¢. 
Since the radius a of the die aperture is given and 


R = a/sin 6, (1) 


the considered stage of the deformation is conveniently specified by the parameter @. 
All “rates’”’ used in the following are rates of change with respect to 6. 

Denote by A the uniform sheet thickness in the bulged state and by h, that in the 
undeformed state. Since the sheet material is assumed to be incompressible, it must 
occupy equal volumina in both states. It follows from this condition that 


h = ho cos” 2 ‘ (2) 
The “‘strain rate” ¢, is thus given by 
1 dh 6 
& =F a6 —tan 5 ; (3) 


™W. T. Koiter, C. B. Biezeno Anniversary Volume, H. Stam, Haarlem, Holland, 1953, p. 232. 
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Since the height of the bulge is 


, (4) 


NIisd 


H = R(1 — cos 6) = atan 


we have 


«, = —H/a. (5) 


Let r be the radius of the parallel circle through A and r, the radius of the correspond- 
ing circle in the undeformed sheet. When the condition of incompressibility is applied, 
not to the entire sheet but only to the portion which is bounded by the circle of radius 
ry) in the undeformed state and by the circle of radius r in the deformed state, the following 


relation is obtained: 








he)" g cos y/2 
ilies ro} , 8 5 = "0 cos 6/2 * (6) 
Now, 

r=Rsing =a~~. (7) 


Eliminating r between (6) and (7) and solving for 7, we obtain 


sin ¢/2 
Tf =a. 8 
, sin 6/2 (8) 





For a given particle ry, remains constant during the deformation. Differentiating (8) 
with respect to @ and using (6) and (8), we therefore obtain 


dp tang/2 7 
— = —=—, 9 
dé tan6/2 = ar (9) 





From (6) and (9) it is seen that the circumferential strain rate is given by 


1 dr d 1 6  o« e| 
= =— = - = — = — = tan” = |. 10 

€. aa 79 (108 r) 5 [tan 2 cot 5 t 5 (10) 

From (3), (10), and the condition of incompressibility, the meridional strain rate is 


then obtained as 
én = —€. — & = [tan 8 + cot 8 tan’ e|. (11) 


The circumferential strain rate vanishes at the edge (g = +6); for relevant values of 
6, say 0 < @ < 2/2, and all values of g between —8@ and 8, the strain rates (10) and (11) 
are seen to be positive. The state of -flow considered here therefore is of the type c), 
and ¢. = dq = ¢. 

4. Static considerations. If o. = o, = ¢, the equilibrium of the assumed spherical 
bulge of constant wall thickness under the applied pressure p requires that o has the 


constant value 





a « Be 12 
a 2he sin 6 cos’ 6/2 (12) 
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in the entire bulge. This means that, at any given stage of the bulging, the entire material 
is in the same state of hardening. 
Consider first a perfectly plastic material that flows under the constant stress o» in 


simple tension. We then have ¢ = o, = const., and hence, from (12), 
fo,hyo . @ , 0 : 
p = —— sin 5 cos’ 5 (13) 


As @ increases, the pressure (13) reaches the maximum value 


ae/ 
3V3 ooh 
max p (14) 
ta 
when 6 60°. The considered spherical bulge of constant wall thickness is not stable 
beyond this pressure maximum. When the pressure maximum is reached, h/ho 3/4 


according to (2). The logarithmic strain in the direction normal to the sheet therefore 
is nearly —30%. 
Consider next a material that strain-hardens according to 


o = 0,(1 +a & }) (15) 


in simple tension, & being the logarithmic strain, and o) and @ being constants. The 


considered state of stress in the bulged sheet has the principal values o, 0 0, 
o, = 0; it may be obtained by the superposition of the state of balanced triaxial tension 
GC. = 6, =G, o and the state of simple compression g. = o, = 0, 6, = —o. The first 


of these states of stress will not produce plastic flow or strain-hardening of the incompres- 
sible material, the second may be assumed to produce strain hardening in accordance 
with (15) prov ided that & in this formula is replaced by log (h/ho), where h/hy is given 
by (2). Substituting 


\ 


0 é 
go=o0 (1 — 2a log cos 5) 16) 


c 





2 Pmox/(ho%), H/o» | 
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into (12) and solving for p, we obtain 


‘ ol 0 . © 
p= toot (1 — 2a log cos é) sin Mu cos* : (17) 


a 


Figure 2 shows the pressure maximum computed from (17) and the corresponding 
values of 0, h/hy , and | &| = —log (h/h»), all versus the strain-hardening parameter a. 
[t is seen that the pressure maximum p,,, as well as the corresponding values of 6 
and | € | increase with a, whereas the ratio h/ho at the pressure. maximum decreases 
with a. In the considered range of @ all these quantities vary with a@ in a nearly linear 


manner 


BOOK REVIEWS 


Mathematical methods for scientists and engineers. By Lloyd P. Smith. Prentice-Hall, Inc., 
New York, 1953. x + 453 pp. $10.00. 


The material in this book is essentially that taught for a number of years by the author to graduate 
students in physics and physical chemistry at Cornell University. The unusual feature of this text is 
the wide range of mathematical methods treated. At the same time it will be found that the treatment 
of the material is quite adequate for dealing with most physical problems. The discussion is consise 
and clear. 

There is no specific treatment of differential equations except as they arise in the treatment of the 
othe I topics 

The reviewer believes that this is one of the most useful books available on intermediate and ad- 
vanced mathematical methods. 

The chapter headings of this text are as follows: Elements of Function Theory; Differential Calculus, 
Integral Calculus; Space Geometry; Line, Surface, and Multiple Integrals; Theory of Functions of a 
Complex Variable Residues and Complex Integration; Representation of Functions by Infinite Series 
of Functions; Applications of Functions of a Complex Variable to Potential and Flow Problems; Algebra 
of Linear Equations, Transformations and Quadratic Forms; Vector and Tensor Analysis; Orthonormal 
Function Systems; Orthonormal Functions with a Continuous Spectrum; Integral Equations; Varia- 
tional Methods; and Elements of Probability Theory. 

Roun TRUELL 


Calculus of variations with applications to physics and engineering. By Robert Weinstock. 
MeGraw-Hill Book Company, Inc., New York, 1952. x + 326 pp. $6.50. 


According to the preface, this volume presents an introduction to the calculus of variations followed 
by application of the subject to problems of physics and theoretical engineering. 

The first five chapters give the usual elementary treatment of the calculus of variations with no 
pretense of complete mathematical rigor. Chapters 6-12 present applications to dynamics, elasticity, 
quantum mechanics, and electrostatics. These applications are, for the greater part, of an elementary 
nature; modern problems in acoustics, electromagnetic theory, and quantam mechanics are not dis- 
cussed 

Up to the present time there is no other volume in the English language that offers such a variety 
of applications of the calculus of variations to problems in physics. The book must therefore be accepted 
as a worthwhile contribution to the applied mathematician’s library. 

This reviewer has a few adverse comments to make on the contents of the book; these are not in- 
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tended to detract from the general acceptability of this volume as a text-book at the first-year graduate 
level. 
(a) It is claimed in the preface that ‘the reader who has mastered the essence of the material included 
should have little difficulty in applying the calculus of variations to most of the subjects which have 
been squeezed out.” This is unwarranted optimism on the part of Professor Weinstock. The book does 
not offer any examples of variational principles for integral equations; there is no mention of Green’s 
functions. How, then, could a reader be expected to understand “with little difficulty” the Levine- 
Schwinger treatment of scattering problems? 

(b) There is no mention of variational techniques for obtaining lower bounds to eigenvalues (work 
of Weinstein, Aronszajn, etc.). In particular the chapter on quantum mechanics is weak because of 
this and other omissions 

(c) The work of Pélya and Szegé and their colleagues on ‘isoperimetric inequalities in mathematical 
physics’’ is barely touched on in the last chapter on electrostatics. 

(d) Chapter 5 is entitled ‘Geometrical Optics: Fermat’s Principle.” This chapter is five pages long 
and consists exclusively of deriving the two-dimensional principle of Fermat by two different methods. 
It is unfortunate that the author did not discuss “corner conditions” in this connection. Snell’s law 
could then be taken as an illustration and the chapter expanded to reasonable size. (See e.g. H. Lewy: 
“Aspects of the Calculus of Variations.’’ University of California Press, Berkeley, California, 1939. 
This excellent little volume is not listed in Weinstock’s bibliography. 

(e) A few serious errors in the text have been pointed out by Synge in his review in ‘‘Bulletin of 
the American Mathematical Society’’, Vol. 59, No. 4, July, 1953. 

Here are two additional non-trivial errors. 

1. The second part of theorem (d), p. 157, is incorrect. 

2. Introduction, p. 2: “The only values of z at which y = g(x) can possibly achieve a minimum 
0.” A minimum ean occur at points where g’(x) does not exist, or at an endpoint 


are the roots of g’(r) = 
of an interval. This and similar precautions are neglected throughout the book; for instance, on p. 39, 


where no condition on the differentiability of g(x, y) is given. 
I. STaKGOLD 


Numerische Behandlung von Differentialgleichungen. By UL. Collatz. Springer-Verlag, 
Berlin, Géttingen, Heidelberg, 1951. xiii + 458 pp. $11.43. 


This book represents a careful and comprehensive study of the numerical methods of solution of 
ordinary and partial differential equations. It also contains a concise treatment of numerical methods 
for integral equations. By virtue of the wealth of material covered the book should prove to be a very 
valuable source of reference to anyone interested in this important subject. 

In addition to a description of the various methods the author has included, particularly in the case 
of ordinary differential equations, analyses of their theoretical foundations. In many instances there is 
a discussion of the estimation of the error. There are numerous illustrative examples, many of them 
representing practical physical problems. 

Chapter I contains, in addition to introductory remarks, difference formulae and quadrature formu- 
lae, a study of initial-value problems for ordinary equations of first and higher order. Chapter IT is 
devoted to boundary-value problems (boundary conditions prescribed at more than one point) for 
ordinary equations, including eigenvalue problems. In Chapter III we find initial-value problems and 
combined initial and boundary-value problems for partial differential equations. In chapter IV the 
author takes up boundary-value problems for partial differential equations. Chapter V is concerned 


with integral equations and a brief discussion of functional equations. 
G. W. Morcan 
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Die praktische Behandlung von Integralgleichungen. By H. Biickner. Springer-Verlag, 
Berlin-Goettingen-Heidelberg, 1952. vi + 127 pp. DM 18.60. 


As the author states in the Introduction, the numerical treatment of integral equations is a compara- 
tively new and still developing field of numerical analysis. The present systematic survey is therefore 
particularly welcome even though it is restricted to equations of the Fredholm type for functions of a 
single variable. The first chapter recapitulates the principal results of the general theory. Chapter II 
presents formulas and variational principles that are useful in the calculation of characteristic values. 
Particular attention is given to methods of bracketing characteristic values. Chapter III is devoted to 
iterative methods. The use of mixed iteration for the solution of inhomogeneous integral equations is of 
special interest. Chapter IV is concerned with techniques using approximations to the kernel of the 
integral equation. Perturbation methods, the methods of Ritz and Galerkin, and methods replacing the 
integral by a sum are treated as special cases. The brief final chapter deals with special] kernels. 


W. PRAGER 


Mathematical physics. By Donald H. Menzel. Prentice-Hall, Inc., New York, 1953 
v + 412 pp. $8.50. 


This text contains five chapters entitled: I. Physical Dimensions and Fundamental Units, II. Me- 
chanics and Dynamics, ITI. Waves and Vibrations, IV. Classical Electromagnetic Theory, V. Relativity. 
The book is claimed by the author to be designed for use in junior, senior, or graduate courses in mathe- 
matical physics; this purpose seems to have been achieved very nicely. The book appears to the reviewer 
to be a very satisfactory text for the transition from advanced undergraduate physics to first year gradu- 
ate physics, and it should be welcomed by many students and teachers. 

This is not a text on mathematical methods in physics. It is rather a text on theoretical physics with 


emphasis on the mathematics. 
Roun TRUELL 


The principles of the control and stability of aircraft. By W. J. Duncan. Cambridge Univer- 
sity Press, 1952. xvi + 384 pp. $8.00. 


In view of the importance of the field and the abundance of texts on other aspects of aeronautics, 
the scarcity of textbooks on aircraft stability and control is surprising. The present volume is therefore a 
particularly welcome addition to the aeronautical literature. Within the limited space available for this 
review, the scope of the work is best indicated by the following chapter headings: Introductory survey. 
Elementary mechanics of flight. The equations of motion of rigid aircraft. Methods for solving the 
dynamical equations and for investigating stability. Longitudinal-symmetric motion. Lateral-antisym- 
metric motion. Flap controls in general. The measurement of aerodynamic derivatives. Controls for 
roll, pitch, and yaw. Static stability and maneuverability. Stalling and the spin (by A. D. Young). 
The influence of distortion of the structure. The influence of the compressibility of the air. Flaps for 
landing and take-off (by A. D. Young). Sundry topics. 

The exposition is particularly clear and easy to follow. 

W. PRaGER 


CORRECTION 


In the review of Petrovskij, Vorlesungen ueber die Theorie Integralgleichungen, this Quarterly 11, 375 
(1953) the price of this book was erroneously stated to be $7.80 instead of DM 7.80. 
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Complex analysis. An introduction to the theory of analytic functions of one complex 
variable. By Lars V. Ahlfors. McGraw-Hill Book Co., New York, Toronto, London, 
1953. xi + 247 pp. $5.00. 


This book is a rigorous introductory treatment of the theory of analytic functions of one complex 
variable. It is very readable and would be an excellent basis on which to teach a course. There are exer- 
cises every few pages. Although no physical examples are given, the student who has mastered this 
book will have no mathematical difficulty applying the knowledge he has gained to any application of 
complex variable in mathematical physics. 

In Chapter 1, complex numbers are introduced, their geometrical representation and linear trans- 
formations considered. Chapter 2 is devoted in the main to theorems needed from the theory of a real 
variable and from the theory of sets. Complex integration is introduced in Chapter 3. Cauchy’s integral 
formula and the calculus of residues follow. Chapter 4 contains convergence and uniform convergence 
of series, Taylor and Laurent series, infinite products and normal families. The Dirichlet problem is 
considered in Chapter 5 and multiple-valued functions in Chapter 6. 


D. R. BLAND 


Proceedings of Symposia in Applied Mathematics—Fluid Dynamics. Volume IV. McGraw- 
Hill Book Co., Ine., New York, Toronto, London, 1953, for the American Math- 
ematical Society, 80 Waterman Street, Providence, Rhode Island. v + 186 pp. $7.00. 


This volume contains thirteen papers and one abstract of a paper which were presented at the 
Fourth Symposium in Applied Mathematics of the American Mathematical Society held at the Univer- 
sity of Maryland on June 22 and 23, 1951. All papers are concerned with fluid dynamics. Topics in 
turbulence, compressible (including transonic), potential, and viscous flow are treated, and one paper 
deals with the development of the hydrodynamical equations from the point of view of thermodynamics 
of irreversible processes. 

The papers included are: 8. Chandresekhar, “Some aspects of the statistical theory of turbulence’. 
C. C. Lin, ‘‘A critical discussion of similarity concepts in isotropic turbulence’. A. Busemann, ‘The 
nonexistence of transonic potential flow’. R. E. Meyer, “(On waves of finite amplitude in ducts (ab- 
stract)’’. T. Y. Thomas, ‘‘On the problem of separation of supersonic flow from curved profiles’’. G. F. 
Carrier and K. T. Yen, “‘On the construction of high-speed flows’. M. H. Martin and W. R. Thickstun, 
‘‘An example of transonic flow for the Tricomi gas’’. A. E. Heins, ‘On gravity waves”. S. R. De Groot, 
“Hydrodynamics and Thermodynamics”. J. M. Burgers, ‘‘Non-uniform propagation of plane shock 
waves”. T. Theodorsen, ‘“‘Theory of propellers’. G. Birkhoff, D. M. Young and E. H. Zarantonello, 
“Numerical methods in conformal mapping’’. J. L. Synge, “Flow of viscous liquid through pipes and 
channels’. A. Weinstein, ‘“The method of singularities in the physical and in the hodograph plane”’. 


GEORGE W. MorGAN 


Anfangswertprobleme bei partiellen Differentialgleichungen. By R. Sauer. Springer- 
Verlag, Berlin-Goettingen-Heidelberg, 1952. xiv + 229 pp. $6.90. 


The book treats initial value problems for single partial differential equations and systems of partial 
differentia] equations of hyperbolic type. The introductory first chapter discusses the typical boundary 
value problem for the Laplace equation and the typical initial value problem for the wave equation. 
The second chapter develops the theory of characteristics for first partial differentia] equations of the 
first order. Chapters III and IV which constitute about three quarters of the volume are devoted to 
the system of quasilinear differential equations of the first order and the equation of the second order, 
Chapter III treating the case of two and Chapter IV the case of more than two independent variables. 

The clear exposition makes the book easy to read. A particularly valuable feature is the constant 
detailed reference to relevant physical problems (hydrodynamics, acoustics, gas dynamics, plasticity). 

W. PRAGER 








A FIRST COURSE IN ORDINARY DIFFERENTIAL EQUATIONS 


By RUDOLPH E. LANGER, Professor of Mathematics, University of Wisconsin. 
In this logical introduction to an important branch of modern mathematics, Dr. Langer 
treats his subject both from the point of view of its mathematical interest and its utility as 
a tool in science and technology. His book stresses a commonsense approach that plays 
down mere memorization and provides a solid understanding of the subject. The author 
emphasizes differential equations of the first and second orders. A useful feature of the 
work is its many problems which permit the reader to develop some familiarity with the 
equations that come up for discussion before he is called on to examine theoretical con- 
siderations. The logical arrangement of the material, which makes it possible to omit 
whole sections without losing continuity, gives the book a lasting value as a reference in 
addition to its immediate use as an introduction to the field. 1954. 249 pages. $4.50. 


THEORY oF GAMES AND STATISTICAL FUNCTIONS 


By DAVID BLACKWELL, Professor of Mathematics, Howard University, and 
M. A. GIRSHICK, Professor of Statistics, Stanford University. Here is a unified method 
for developing statistical decision concepts from the point of view of game theory. It is the 
author’s thesis that that decision theory applies to statistical problems the principle that a 
statistical procedure should be evaluated by its consequences in various circumstances. The 
mathematical model for a decision theory is a special case of that for game theory, and in 
this book relevant parts of game theory are used as bases for a self-contained, rigorous ex- 
position of statistical decision theory. The book offers a new, unified method of developing 
statistical concepts affording a clearer insight into the problems of design and analysis of 
experiments. A wide selection of problems is included. 1954. 356 pages. Probably $7.00. 


THEORY oF EQUATIONS 


By CYRUS COLTON MACDUFFEE, Professor of Mathematics, University of Wis- 
consin. This multi-purpose text is planned to meet the needs of two distinct groups of 
students: those who will go on from the basic course to graduate work in mathematics, 
and those who seek only the mathematical skills they require as “tools” in pursuing their 
own particular branches of science. The work meets these divergent requirements by cover- 
ing the standard material in a fairly conservative manner, and at the same time introducing 
the important concepts of modern algebra which enable the graduate to go on to abstract 
algebra without dislocation. The book gives a complete treatment of linear systems in terms 
of the modern methods which avoid the use of determinants. The book also provides an 
unusually thorough discussion of systems of equations of higher degree, and a concise in- 
troduction to the theory of numbers. 1954. 120 pages. $3.75. 
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